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1  Introduction 


An  effort  aimed  at  improving  the  quality  of  the  numerical  analysis  of  complex  flow-thermal- 
structural  interactions  in  hypersonic  Bight  was  motivated  by  a  recent  surge  of  interest  in 
developing  new  supersonic  and  hypersonic  vehicles  (as  in  the  National  Aerospace  Plane 
(NASP)).  Such  structures,  in  particular  their  leading  edges  and  points  of  shock  impingement, 
are  exposed  to  hostile  aero-thermomechanical  loads  with  resulting  elevated  temperatures. 
Even  when  equipped  with  active  cooling  systems,  they  can  reach  temperatures  close  to  the 
melting  point.  For  a  designer,  these  structures  present  a  significant  challenge  because  of  the 
inherent  nonlinearities  in  all  aspects  of  the  multidisiplinarv  analysis.  Some  of  the  critical 
computational  issues,  identified  in  Ref.  [S],  include  the  difficulties  associated  with:  (1) 
analyzing  the  viscous,  compressible  flow  and  predicting  the  high  local  aerodynamic  heating. 
(2)  modeling  and  analyzing  of  multimode  unsteady  heat  transfer  in  a  high  temperature 
convectively-cooled  structure,  and  (3)  simulating  the  transient,  nonlinear  thermal-structural 
response  to  rapid  temperature  changes. 

The  general  objective  of  this  project  was  to  make  further  progress  toward  a  realistic 
analysis  and  design  of  hypersonic  structures.  The  effort  was  threefold: 

•  identification  of  the  major  physical  phenomena  involved  in  flow-thermal-structural  in¬ 
teractions  in  hypersonic  flight. 

•  development  of  computational  methods  and  techniques  necessary  for  reliable  modeling 
and  analysis  of  these  phenomena,  and 

•  solution  of  examples  of  hypersonic  structures  subjected  to  severe  aerothermal  loads 
due  to  fluid-thermal-structural  interactions. 

In  the  course  of  this  project  significant  progress  has  been  made  in  all  these  areas  and 
the  results  are  presented  in  this  report.  A  general  summary  of  the  research  progress  is 
summarized  in  Section  2.  A  theoretical  formulation  and  finite  element  approach  for  the 
solution  of  transient  thermal  problem  is  presented  in  Section  3.  This  section  is  followed 
by  a  finite  element  formulation  of  tne  thermo-viscoplastic  structural  analysis  with  large 
displacements  (Section  4).  In  Section  5  a  derivation  and  finite  element  formulation  of  a 
new  implicit/explicit  Taylor-Galerkin  method  for  the  solution  of  compressible  flow  and  the 
prediction  of  aerothermal  loads  is  presented.  In  Section  6  several  validation  and  test  examples 
for  thermo-structural  and  flow  analyses  are  solved.  Finally,  in  Section  7,  a  more  realistic 
thermo-viscoplastic  analysis  of  a  segment  of  the  scramjet  strut  is  presented.  This  section  is 
followed  bv  conclusions  and  final  remarks. 
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2  Summary  of  Research  Progress 


The  work  in  the  project  has  focused  on  the  physical  phenomena,  mathematical  formulations 
and  computational  methods  needed  to  solve  fluid-thermal-structural  interactions  in  hyper¬ 
sonic  flight.  Toward  this  end,  unique  computational  capabilities  and  finite  element  software 
have  been  developed. 

The  flow-thermal-structure  (FTS)  interaction  problem  considered  here  is  shown  in  Figs.  1 
and  2.  Figure  2a  shows  a  low  temperature,  undeformed  elastic  structure  subjected  to  a  low- 
level.  uniform,  aerodynamic  heating.  In  this  state,  there  is  little  FTS  interaction.  In  contrast. 
Fig.  2b  shows  a  structure  that  is  subjected  to  intense  localized  heating  characteristic  of  shock 
interactions.  There  is  a  local,  growing  region  of  very  high  temperatures  and  plastically 
deforming  material.  As  the  plastic  region  grows,  the  surface  deforms  into  the  flow  producing 
flow  recirculation.  The  flow  recirculation  intensifies  local  heating,  skin  friction,  and  surface 
pressure.  As  the  heating  increases,  the  temperature  continues  to  grow  and  can  eventually 
reach  the  melting  temperature. 

Toward  a  realistic  analysis  of  the  above  FTS  interaction  problem  the  following  objectives 
have  been  accomplished: 

1.  Basic  phenomena  essential  in  the  analysis  of  these  problems  have  been  identified.  For 
these  phenomena,  mathematical  models  have  been  formulated,  assuming  that: 

(a)  The  unsteady  heat  transfer  problem  is  described  by  conservation  of  energy  ne¬ 
glecting  mechanical  coupling  terms.  This  assumption  is  known  to  be  reasonable 
for  the  materials,  deformation  rates  and  applied  heat  loads  encountered  in  hyper¬ 
sonic  structures.  Nonlinear  heat  transfer  due  to  temperature  dependent  material 
properties  such  as  specific  heat  and  thermal  conductivity  are  included  as  well  as 
surface  radiation. 

(b)  The  structural  problem  is  assumed  to  be  quasi-static  and  viscoplastic  with  large- 
deformation  kinematics.  The  application  of  a  viscoplastic  constitutive  model  is 
necessary  because  of  the  extremely  high  temperatures  enocuntered  in  hypersonic 
structures.  The  particular  constitutive  theory  applied  is  known  as  the  Bodner- 
Partom  unified  viscoplastic  theory  [4 1 .  This  theory  has  been  validated  lor  high- 
temperature  alloys  under  the  NASA  HOST  program  [5.6].  An  extension  oi  this 
theory  for  nonlinear,  large  displacement  kinematics  was  necessary  in  order  to 
capture  beam-column  effects  occurring  due  to  localized  heating  (see  Section  7). 

(c)  The  external  flow  problem  is  governed  by  the  Navier-Stokes  equations  for  a  com¬ 
pressible.  viscous  flow  of  a  calorically  perfect  gas.  This  assumption  is  legitimate 
for  supersonic  and  hypersonic  flows  where  the  Mach  number  is  less  than  approx- 


imately  eight.  For  higher  Mach  numbers,  when  the  gas  temperature  is  2500  K 
or  higher,  the  specific  heat  becomes  temperature  and  pressure  dependent  and  the 
air  becomes  chemically  reactive. 

Computational  methods  and  finite  element  software  have  been  developed  for  the  an¬ 
alysis  of  flow-thermal-structural  interaction  problems  described  by  appropriate  math¬ 
ematical  models,  in  particular: 

(a)  A  family  of  finite  element  heat  transfer  programs  has  been  placed  into  operation. 
These  programs  solve  nonlinear,  steady  or  transient  heat  transfer  problems,  taking 
into  account: 

•  conduction  with  temperature-dependent  properties 

•  surface  convection 

•  convective  cooling  (active  cooling  systems) 

•  radiation 

(b)  A  two-dimensional  finite  element  thermo-viscoplastic  structural  analysis  program 
has  been  developed.  The  program  is  coupled  with  the  thermal  analysis  software 
and  performs  a  nonlinear,  transient  solution  of  structural  deformation  problems, 
taking  into  account: 

•  visco-elastoplastic  material  models  with  temperature-dependent  properties 

•  large  displacement  kinematics  (to  capture  beam-column  effects  due  to  local¬ 
ized  heating) 

•  development  of  continuum  damage  in  the  material  (to  predict  structural  fail¬ 
ure) 

The  above  viscoplastic  problems  result  in  a  set  of  nonlinear.  "  stiff"  evolution 
equations.  In  order  to  provide  a  reliable,  efficient  solution  for  these  problems, 
the  finite  element  software  was  equipped  with  an  error  estimator  and  an  adaptive 
time  stepping  procedure  which  enables  precise  control  of  the  solution  error. 

(c)  A  new,  implicit/explicit  method  for  the  compressible  Xavier-Stokes  equations  was 
formulated  and  implemented  in  the  adaptive  finite  element  program.  The  method 
was  developed  for  the  prediction  of  aerothermal  loads  on  hypersonic  structures. 
I’hc  development  of  a  new  implicit/explicit  method  was  motivated  by  our  numer¬ 
ical  experiments  and  the  experience  of  other  researchers  that: 

•  Explicit  solvers  are  inefficient  in  the  solution  of  viscous  supersonic  (lows  due 
to  time  step  stability  limitations  and  slow  convergence  of  the  heating  rates. 

•  Fully  implicit  solvers,  while  being  more  stable  and  reliable  than  explicit  al¬ 
gorithms.  require  large  computational  times  due  to  expensive  element  matrix 
evaluations  and  solution  of  linear  system  of  equations. 
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The  new  implicit/explicit  method  developed  in  the  project  adaptively  selects  im¬ 
plicit  and  explicit  computational  subdomains  to  provide  stability  of  the  implicit 
methods  while  at  the  same  time  minimizing  the  cost  of  the  computations. 

(d)  Numerical  solutions  of  a  variety  of  problems  have  been  performed.  These  problems 
include: 

•  relatively  simple  test  problems,  designed  to  verify  methods  and  procedures 
and  compare  results  with  existing  experimental  evidence  and  analytical  solu¬ 
tions.  and 

•  more  realistic  interaction  problems,  simulating  the  behavior  of  aerospace 
structures  under  severe  aerothermal  loads  in  hypersonic  flight. 

It  should  be  noted  that  the  major  effort  in  the  first  two  years  of  the  project  was  oriented 
toward  identification  ot  physical  phenomena  in  hypersonic  interactions,  development 
of  corresponding  mathematical  models,  and  computational  procedures  and  their  nu¬ 
merical  verification.  Thus,  the  number  of  realistic  examples  solved  in  this  phase  is 
rather  small.  More  extensive  computations,  combined  with  experimental  verification, 
was  planned  for  the  third  year  of  the  project. 

Papers  and  Presentations 

The  research  performed  in  the  project  and  the  corresponding  results  were  presented  at  the 
AIAA/ASME/ASCE/AMS  30th  Structures,  Structural  Dynamics  and  Materials  Conference 
in  Mobile.  Alabama  in  April  19S9  (Ref.  [33]).  The  corresponding  paper  was  accepted  for 
publication  in  the  AIAA  Journal  in  1989.  Ref.  [35].  Further  papers  and  presentations  are  in 
preparation. 


3  Thermal  Analysis  of  Hypersonic  Structures 

Convert  ively  cooled  structures  are  strong  candidates  for  use  as  hypersonic  flight  vehicles.  For 
hypersonic  flight,  some  leading  edges  and  panels  require  active  cooling  systems  to  keep  the 
structural  temperature  within  acceptable  ranges.  The  internal  flow  in  the  coolant  passage 
pla_  ,  a  predominant  role  in  the  thermal  response  of  a  hypersonic  structure  subject  U>  externa! 
heating.  A  cross-section  of  a  typical  convectively  cooled  structure  is  shown  in  Fig.  3.  An 
aerodynamic  skin  and  coolant  passage  with  an  internal  heat  exchanger  protect  the  primary 
structure  from  the  aerodynamic  heating.  The  thin,  typically  metallic,  aerodynamic  skin 
transfers  the  energy  from  the  aerodynamic  heating  to  a  low  temperature  coolant  flow  through 
the  heat  exchanger  fins  that  connect  the  aerodynamic  skin  to  the  primary  structure.  In  a 
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typical  engine  structure,  the  coolant  is  cold  hydrogen  that  later  is  used  as  the  propulsion 
system  fuel. 

Heat  transfer  in  the  aerodynamic  skin  consists  of  conduction  combined  with  surface 
radiation.  Heat  transfer  between  the  aerodynamic  skin,  the  heat  exchanger  surfaces  and  the 
primary  structure  is  by  conduction  at  the  solid-fluid  interface.  The  dominant  mode  of  heat 
transfer  in  the  coolant  flow  is  forced  convection. 

An  accurate  representation  of  the  heat  transfer  in  the  coolant  passage  is  the  critical  step 
in  the  heat  transfer  analysis.  There  are  two  basic  representations  that  can  be  employed.  The 
first,  denoted  here  as  the  engineering  model,  is  based  upon  a  number  of  assumptions  that 
greatly  simplify  the  problem  into  a  single  energy  equation  with  a  specified  mass  flow  rate. 
Detailed  computations  of  the  fluid  velocity  components  and  temperature  is  not  required. 
The  second  representation  idealizes  the  coolant  flow  as  a  continuum  model,  and  the  partial 
differential  equations  describing  conservation  of  mass,  momentum  and  energy  are  solved 
simultaneously  to  obtain  fluid  velocity  and  temperature  distributions.  This  latter  model  is 
the  most  accurate,  but  it  is  also  considerably  more  expensive  than  the  engineering  model. 
As  a  first  step  towards  understanding  flow-thermal-structural  interactions  characteristic  of 
hypersonic  structures  the  engineering  heat  transfer  model  is  employed. 

3.1  Engineering  Model  for  Coolant  Passages 

The  basic  features  of  the  engineering  model  can  be  developed  from  the  idealization  shown 
in  Fig.  4.  A  segment  of  the  coolant  passage  of  width  w  is  shown:  for  simplicity,  only  the 
upper  one-half  of  the  coolant  passage  with  the  aerodynamic  skin  is  shown.  The  engineering 
formulation  (Ref.  [34])  is  based  on  the  following  assumptions: 

1.  The  thermal  energy  state  of  the  fluid  is  characterized  by  the  fluid  bulk  temperature 
Tp  which  varies  only  in  the  flow  direction,  i.e.  Tp(x.t). 

2.  The  flow  is  represented  by  the  mass  flow  rate  m  in  the  coolant  passage  specified  by 
in  =  ppApVp  where  pp  is  the  coolant  density,  Ap  is  the  cross-sectional  area  ot  the 
coolant  passage,  and  Vp  is  the  coolant  mean  flow  velocity. 

3.  A  convection  coefficient  h  is  defined  such  that,  the  heat  flux  qp  transferred  between  the 
structure  and  t  he  coolant  may  be  expressed  as 

q  =  h(Ts  -  Tp) 

where  T^(x.t)  denotes  the  structural  temperature  at  the  fluid-solid  interface. 
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4.  The  convection  coefficient  h  may  be  expressed  as  a  function  of  the  fluid  bulk  temper¬ 
ature  alone  by  using  the  analytical/empirical  equation  for  the  Mussel t  number. 

_ ,  hD 

-V“  =  I7 

where  D  is  the  hydraulic  diameter  of  the  coolant  passage,  and  kF  is  the  thermal 
conductivity  of  the  fluid  coolant.  The  convection  coefficient  as  well  as  the  solid  and 
fluid  thermal  parameters  may,  in  general,  be  temperature  dependent. 


W  ith  these  assumptions,  energy  balances  on  the  aerodynamic  skin  and  coolant  give  the 
governing  conservation  equations. 

Fluid: 
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Solid: 

c)  d'T  IT 

-  fc{ksAs~Q^r)  +  wfl(Ts  -  Tf)  +  creT4  +  pscsAs—jA  =  q  (3.2) 

where  the  subscripts  F  and  5  denote  the  fluid  and  solid,  respectively.  In  these  equations. 
c  denotes  specific  heat,  cr  is  the  Stefan-Boltzmann  constant,  and  e  is  the  surface  emissivity. 
Because  ot  the  temperature  dependence  of  the  thermal  parameters  and  the  radiation  term. 
Eqs.  (3.1)  and  (3.2)  constitute  a  nonlinear  set  of  partial  differential  equations. 

The  heat  exchanger  fins  are  not  included  explicitly  in  this  model.  However,  the  heat 
transfer  between  the  heat  exchanger  fins  can  approximately  be  taken  into  consideration 
though  the  use  of  an  effective  width  w  which  represents  the  area  over  which  the  convective 
heat  exchange  occurs. 


3.2  Finite  Element  Formulation 

A  typical  finite  element  representing  Eqs.  (3.1)  and  (3.2)  is  characterized  by  fluid  and  tluid- 
solic!  interface  nodes.  The  element  shown  in  Fig.  5  has  two  fluid  nodes  (I  and  .J).  and  two 
fluid-solid  interlace  nodes.  Within  the  element,  the  fluid  and  solid  temperatures  are  express 
as 

Tf  =  [-V(x)]{7>} 


Ts  =  [X(x)){Ts}  (3.3) 

where  f.V(.r)]  are  the  element  interpolation  functions.  Following  usual  finite  element  proce¬ 
dures  (Ref.  [34]),  the  discritized  equations  for  an  element  of  length  L  may  be  derived  in  the 
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form 


- , 

o 

■’1 

o 

Tf 

A  „  +  A  *  +  A  f  —  A 

1 

.  + 

0  Cs 

_  _j 

1  TS  J 

-  AV  AV  +  AV  . 

where  the  element  capacitance  matrices  are  given  by 

[C/r]  =  f  PfCfAp{N}[N}(1x 

Jo 


[Cs]  =  [L  pscsAs{N}[N}dx 
Jo 

matrices  are 

rL  dN 

=  J  mcF{N}[—}dx 


and  the  element  conductance 

[A',1 

[AV] 

[AVI 

[AV] 


Q  J 


(3.5) 


(3.6) 


[I<r]{Ts}  =  j"atTA{S}dx 

The  element  equations  given  in  Eqs.  (3.4)  show-  that  the  coolant  passage  model  can  be 
regarded  as  an  assembly  of  elements  where  each  element  represents  a  single  heat  transfer 
mode.  Thus  we  can  represent  the  coolant  passage  in  terms  of  two  convective  elements:  (T )  a 
mass-transport  element  (Fig.  6a).  and  (2)  a  surface  convection  element  (Fig.  6b).  1  he  mass 
transport  element  represents  the  downstream  convective  heat  transfer  due  to  mass  How:  it 
lias  a  capacitance  matrix  [CF\  as  well  as  conductance  matrices  [AV]  and  [AV].  The  surface 
convection  element  represents  the  convection  heat  exchange  between  the  coolant  fluid  and  the 
solid:  it  does  not  have  a  capacitance  matrix,  but  it  has  a  conductance  matrix  [ A V j  linking  the 
solid  and  fluid  nodes.  These  coolant  passage  elements  are  assembled  into  the  finite  element 
thermal  model  for  the  complete  convecavely  cooled  structure  that  has  conduction  elements, 
radiation  elements,  and  surface  elements  for  convection  to  a  specified  convective  exchange 
temperatu  re. 

The  unsteady  thermal  analysis  is  nonlinear  because  of  the  temperature  dependent  thermal 
properties  and  surface  radiation.  The  equations  are  solved  by  time  marching  with  the  (,'rank- 
Xicolson  algorithm:  at  each  time  step,  the  nonlinear  algebraic  equations  are  solved  using 
Newton-Raphson  iteration. 
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4  Thermo— Viscoplastic  Structural  Analysis 


In  this  section  we  will  present  a  new,  extended  formulation  of  the  thermo-viscoplastic  finite 
element  structural  analysis.  Extensions  of  the  theory  compared  with  our  previous  reports 
include: 

•  Generalization  of  the  formulation  for  the  case  of  large  displacements. 

•  Incorporation  of  a  model  of  material  damage  in  the  constitutive  equations. 

The  extension  of  the  theory  to  the  case  of  large  displacements  was  motivated  by  the 
analysis  of  a  convectively  cooled  structure  presented  in  Refs.  [33,35].  In  this  analysis,  which 
simulated  the  failure  of  heat  exchanger  fins,  the  aerodynamic  skin  experienced  relatively 
large  displacements.  Modeling  of  these  displacements  and  the  occurrence  of  a  nonlinear 
beam-column  type  effect  motivated  our  effort  to  extend  the  theory  to  large  deformations. 

The  particular  version  of  the  large  deformation  theory  applied  here  can  be  classified  as 
a  large  displacement-small  strain  theory.  The  limitation  to  small  strains  results  primarily 
from  the  fact  that  the  eiasto-viscoplastic  Bodner-Partom  constitutive  equations  were  nor 
devised  and  validated  for  stretch  levels  exceeding  a  few  percent  (say  five  percent).  Thus  it 
is  reasonable  to  limit  our  theory  and  applications  to  small  strains. 

The  incorporation  of  a  material  damage  model  in  the  constitutive  equations  allows  for  a 
more  complete  analysis  of  the  structure  and  prediction  of  structural  failure  due  to  extensive 
loads  and/or  multiple  load  cycles.  The  model  presented  here  is  based  on  the  work  o(  Rodner 
and  Chan  (Ref.  [4])  and  is  consistent  with  the  Bodner-Partom  model  of  elasto-viscoplasticity. 

4.1  Large  Displacement-Small  Strain  Kinematics 

4.1.1  Strain  and  Stress  Measures 

In  this  section  we  present  the  basic  kinematic  relations  for  the  theory  of  large  deformations 
and  the  implications  of  the  restriction  to  the  small  strain  case. 

We  are  interested  in  the  deformation  of  a  metrical  continuum,  which  at  a  certain  initial 
time  t  =  0  occupies  an  initial  configuration  which  is  formally  identified  with  the  reference 
configuration  Br.  The  reference  configuration  is  parametrized  by  the  reference  (.'artesian 
coordinate  system  { _ V /v- }  with  base  Ik,I\  =  1,2  (see  Fig.  7). 

The  reference  coordinates  X;  identify  particles  of  a  continuum  in  the  sense  that  by 
particle  X  we  mean  the  particle,  which  in  the  reference  configuration  has  a  position  defined 
bv  vector  X. 


S 


At  an  arbitrary  time  t  >  0  the  continuum  occupies  a  deformed  configuration  Bt,  with  new 
positions  of  particles  referred  to  the  Cartesian  coordinate  system  {.rr }  with  base  k  =  1.2. 
It  is  convenient,  without  restricting  the  generality  of  the  description,  to  identify  coordinate 
systems  { x* }  and  {A’/v}  (see  Fig.  7).  The  notion  of  systems  {A\}  and  {x*}  represents 
a  typical  notation  convention  of  large  deformation  kinematics  (see  references  [15.23. 30]  in 
which  uppercase  symbols  and  indices  refer  to  the  reference  configuration,  while  lowercase 
symbols  and  indices  refer  to  the  current  configuration. 

Moreover,  in  general  curvilinear  coordinate  systems,  superscripts  denote  contravariant 
components  and  subscripts  denote  covariant  components  of  vectors  or  tensors.  However, 
in  our  case  both  coordinate  systems  {A'/}  and  {x, }  are  Cartesian,  so  there  is  no  defer¬ 
ence  between  covariant  and  contravariant  coordinates,  and  thus  only  subscripts  will  be  used 
throughout  this  work. 

Motion  of  the  continuum  is  described  by  the  relation 

x,  =  x,(A’/.  0  t  >  0  (-1.1 ) 

and  the  information  about  local  deformation  in  the  neighborhood  of  a  particle  X  is  presented 
by  the  deformation  gradient  tensor  F  with  components: 

F'j  =  Xij 

where  a  comma  denotes  partial  differentation. 

To  indicate  clearly  implica  ions  of  the  restriction  to  small  strain  analysis,  it  is  useful  to 
apply  the  polar  decomposition  theorem  (see  e.g..  references  [15.23.36]).  according  to  which 
the  deformation  gradient  tensor  F  can  be  decomposed  as: 

F  =  RU 

where  U  is  the  right  stretch  tensor  (symmetric,  positive  definite)  and  R  is  a  rotation  tensor. 
The  above  equation  implies  that  the  general  deformation  can  be  viewed  as  a  two-stage 


process 


jr  —  pV-) 


with  tin  first  stage  being  stretch  (generally  arbitrarily  large),  and  the  second  stage  being 
the  rigid  rotation.  Within  our  theory  we  will  limit  the  analysis  to  the  case  when  F{ "  =  U 
represents  infinitesimal  deformation,  while  tire  rigid  rotation  F{~)  =  R  is  arbitrarily  large. 
From  the  above  we  conclude,  by  standard  arguments  of  the  infinitesimal  theory,  that  at  the 
first  stage  the  stress  and  strain  measures  corresponding  to  the  large  deformation  theory  and 
the  infinitesimal  theorv  can  be  identified 


In  the  above,  Eu  and  etJ  are  the  components  of  the  Green  and  Cauchy  strain  tensor,  re¬ 
spectively,  T;j  and  ttJ  are  the  components  of  the  second  Piola-Kirchhoff  and  Cauchy  stress 
tensors,  respectively,  and  , crtJ  are  generic  symbols  of  strain  and  stress  in  the  infinitesimal 
theory. 

Equation  (4.2)  suggests  that  at  the  first  stage  of  deformation  the  constitutive  equations 
can  equivalently  be  applied  to  the  components  Tu  ,E[j  as  well  as  to 

With  regard  to  the  second  stage  of  deformation  F ^  =  R.  it  can  be  easily  verified 
that  it  changes  components  ttJ  and  etJ,  but  components  Tu  and  Eu  are  invariant  under 
rotation  and  they  still  represent,  the  infinitesimal  stretch  of  the  first  stage  =  U  (but 
now  material  fibers  are  rotated). Thus,  for  the  combined  deformation  F  =  u-e  can 

use  the  elasto-viscoplastic  constitutive  equations  expressed  in  terms  of  Tij.Eij.  This  will 
guarantee  satisfaction  of  the  observer  independence  condition,  while  still  allowing  for  the 
application  of  the  small  strain  constitutive  equations. 

Note,  however,  that  for  this  total  deformation  the  components  of  the  strain  tensor  must 
be  expressed  by  the  fully  nonlinear  formula 

Eij  =  +  «j./  +  Uk.l  ukj)  i  =  I.  j  =  J  (4.3) 

One  of  the  consequences  of  the  small  strain  assumption  is  that 

det  U  1  ( 4.4 ) 

which,  by  properties  of  the  rotation  tenser,  gives 

det  F  =  det  R  det  U  1  (4.5) 

4.1.2  Interpretation  of  Strain  and  Stress  Measures 

The  important  issue  in  practical  computations  is  the  physical  interpretaion  of  strain  and 
stress  measures.  Although  the  reference  measures  Tu  and  Eu  have,  in  general,  no  straight¬ 
forward  physical  meaning,  it  can  be  shown  that  in  our  case  they  actually  have  a  clear 
interpretation. 

From  the  equation  (4.2)  valid  for  the  stage  .F*1'  =  U  (infinitesimal)  and  the  invariance 
of  tensor  E  under  rotation  F(2)  =  R  we  can  easily  conclude  that  components  Eu  retain  a 
typical  interpretation  of  the  infinitesimal  theory.  In  particular.  Eu  represents  the  relative 
stretch  of  material  fibers  which  were  originally  parallel  to  axis  A’/  : 


10 


while  Eu ,  with  I  ^  J.  represents  half  of  the  change  of  the  angle  between  such  fibers.  Note 
however,  that  these  fibers  undergo  a  large  rotation  (see  Fig.  7).  To  present  an  interpretation 
of  the  components  of  the  Second  Piola-Kirchhoff  stress  tensor  Tjj.  let  us  introduce  the 
convected  coordinate  system  {(9/},  which  at  the  initial  time  t  =  0  is  identified  with  the 
reference  system  {A'/},  and  at  further  stages  of  motion  deforms  with  the  body,  so  that 
coordinates  of  material  particles  in  this  system  remain  unchanged 

9,(X.t)  =  0,(X,0)  =  X[(X) 


The  convected  coordinate  system  is.  in  general,  curvilinear  and  its  basis  vectors  are  denoted 
by  g{.  I  =  1,2.  The  components  ttJ  of  a  Cauchy  or  "real'’  stress  tensor  in  the  Cartesian 
coordinate  system  {.x,}  are  expressed  via  its  components  tXJ  in  {0;}  by  the  transformation 
law  (since  {0[}  is  curvilinear,  we  will  temporarily  distinguish  covariant  and  contravariant 
components  of  tensors): 


,  in  _  3x'  AJ 

13  dd,  dOj 


dxj  dXj  rj 

dX[  dXj 

where  differentiation  dxJdXj  is  understood  with  reference  to  relation  (4.1). 


(4.6) 


On  the  other  hand,  the  relation  between  Cauchy  and  the  second  Piola-Kirchholf  stress 
tensors  is  by  definition 

=  t'3  =  (detF)~*  ffe'  dx^TIJ  (4.71 

dX,d Xj 

By  comparison  of  formulas  (4.6)  and  (4.7).  and  recalling  equation  (4.5)  we  obtain  the  simple 
relation 

Tu  _  fiJ 


If  we  recall  that  the  system  {A’/}  is  orthonormal  we  can  also  write: 


Tu  =  t1J. 


This  means  that  in  the  case  under  consideration  the  components  of  the  second  Piola-Kirchhoff 
stress  tensor  can  be  interpreted  as  the  components  of  the  Cauchy  stress  tensor  ("real”  stress) 
presented  in  the  base  of  the  convected  coordinate  system  (see  Fig.  7). 


4.2  Elasto— Viscoplastic  Constitutive  Model  With  Damage  and 

Temperature  Effects 

We  now  describe  the  Bodner-Partom  constitutive  equations  used  in  the  thermo-viscoplast  ic 
analysis.  The  version  of  the  equations  presented  describes  nonelastic  behavior  of  mat  ('rials 
including 
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•  viscoplastic  effects, 


•  temperature-dependent  properties, 

•  development  of  internal  damage  in  the  material 

When  combined  with  the  strain  and  stress  measures  introduced  in  the  previous  section,  these 
equations  are  applicable  in  the  large  displacement-small  strain  theory. 

It  is  important  to  note  that  while  we  will  be  using  large  deformation  strain  and  stress 
measures  Eu  and  T[j  in  the  following  developments,  the  form  of  constitutive  equations  is 
exactly  the  same  as  in  infinitesimal  theory,  and  it  is  not  applicable  at  large  stretch  levels. 
Although  one  can  formally  perform  computations  at  large  strain  levels,  the  results  would 
most  probably  be  inadequate  for  modeling  real  physical  phenomena. 

The  elastic-viscoplastic  analysis  with  temperature  effects  is  based  on  decomposition  of 
strain  rates 

Eu  =  E\e)  +  E\n]  (4.S) 

where  superscripts  (e)  and  (n)  denote  elastic  and  nonelastic  strain  components,  respectively. 
The  total  strain  rates  are  expressed  in  terms  of  displacement  rates  as 

E/j  =  -(fq.j  -f  ujj  +  uk.jukj  A  UkjUk.j)  i  =  J-j  —  J 

and  the  constitutive  relations  for  the  elastic  part  are 

Tu  =  CijklEkl  ~  CijklEklT  (-*•-)) 

where  C  is  the  Hooke's  tensor  of  elasticity  parameters.  £  is  a  tensor  of  thermal  expansion 
and  T  represents  the  rate  of  the  change  in  temperature.  Both  C  and  £  are  temperature 
dependent. 

A  nonelastic  deformation  is  governed  by  the  flow  rule: 

E\j]  =  fu(Tu,Zk^k) 

Z,  =  g,(T[j,Zk) 

>*' I  =  h  ,{T[J ,  a.’,  ) 

where  f/j.g,  and  h,  are  constitutive  functions,  Z,  are  internal  state  variables,  and  II,  are 
damage  variables.  These  functions  and  state  variables  characterize  the  viscoplastic  response 
of  the  material  with  continuum  damage  effects. 
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In  the  particular  version  of  the  Bodner-Partom  theory  applied  in  this  work,  the  nonelastic 
how  rule  is  of  the  form: 


r-(n)  _ 


=  AS 


u 


A  >  0 


/?(n)  _  n 

J-JKK  ~  u 

where  Su  are  the  deviatoric  components  of  a  stress  tensor 


Su  =  Tjj  -  -T[j6[j 


The  current  value  of  parameter  A  is  given  by 

1 


=  T  Dq  exp 
■h 


Z2(  1  -S-)\n 
3  J, 


where  J_>  is  the  second  invariant  of  a  deviatoric  stress  tensor 

1 


■h  =  -srJs 


u 

Dq  is  a  limiting  strain  rate  in  shear,  n  is  a  material  constant  and  Z  and  ^  are  state  variables 
which  evolve  during  deformation.  In  particular.  Z  is  the  hardness  variable,  which  represents 
viscoplastic  hardening  (or  softening)  of  a  material.  The  new  variable  not  present  in  the 
previous  versions  of  the  Bodner-Partom  equations,  is  the  damage  variable.  This  variable 
represents  weakening  of  the  material  due  to  nucieation  and  propagation  of  microscopic  voids 
and  cracks  in  the  material.  The  variable  tc  evolves  from  0  to  1.  where  a  value  ^  =  1  represents 
total  loss  of  the  strength  of  the  material. 

In  the  framework  of  materials  science  the  value  of  is  usually  interpreted  as  ratio  ot  the 
area  of  voids  to  the  total  area  of  a  certain  cross  section  of  a  sample: 

_  -Avoid 

^  ~  A 

The  state  variables  Z  and  evolve  according  to  the  specific  equations  ot  the  viscoplastic 
theory: 

1.  Evolution  equations  of  hardness  variable 

The  internal  state  variable  Z  consists  of  isotropic  and  directional  components. 

Z  =  Z7  +  ZD 

The  evolution  equation  proposed  for  the  isotropic  hardening  component  (Reis.  [33. 35]) 


is 


ZI{t)  =  ml[Zx  -  Z/(0]%(<)  -  -AiZ, 


1ri 
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Z!{t)  —  Zi 
Zr 


d.lO) 


with  the  initial  condition.  Z\ 0)  =  Zq.  In  the  first  term,  Z\  is  the  limiting  (saturation)  value 
of  Zl ,mx  is  the  hardening  rate,  and  the  plastic  work  rate  is 

W,  =  TUE\P] 

which  is  taken  as  rhe  measure  of  hardening.  Z2  is  the  minimum  value  of  Zl  at  a  given 
temperature,  and  .  T  and  rl  are  temperature  dependent  material  constants.  The  evolution 
form  of  the  directional  hardening  component  (Refs.  [33.35])  is  defined  as 

ZD(t)  =  3 [j{t)uu(t ) 


where  tz/j  are  the  direction  cosines  of  the  current  stress  state. 

«/j(o  =  Tu^iyrKiTxL)1  (4.1D 

The  evolution  equation  for  3u(t)  has  the  same  general  form  as  that  for  isotropic  hardening 
but  has  tensorial  character. 


3u  =  m2[Z3uu(t)  -  0u(t)}Wp(t) 


—  a2z. 


where 


and 


vu(t)  =  3u{t)/[i3Ki(t)^ix-i(t))i 


3,  A  0)  =  0 


As  in  Eq.  (4.10),  m*  is  the  hardening  rate.  .4 2  and  r2  are  temperature  dependent  material 
constants. 

2.  Evolution  of  damage 

The  damage  parameter  consists,  in  general,  of  isotropic  and  directional  components. 


—  wu  T  U, 


D 


The  evolution  of  isotropic  damage  proposed  in  reference  [4]  is  of  the  form 


In  the  above  p  and  H  are  material  constants  (temperature  dependent),  0  is  the  stress 
intensity  function,  given  by 


Q  = 


2  +  cif 


14 


where  cr,|iax  is  the  maximum  principal  tensile  stress,  /,+  is  the  first  stress  invariant  (nonneg¬ 
ative)  and  J2  is  the  previously  introduced  second  invariant  of  deviatoric  stress. 

,-l.  B.  D ,  and  u  are  material  constants.  A.B.C  must  satisfy  the  condition 

.4  +  B  +  C  =  1 


Clearly,  the  actual  proportion  of  these  constants  selects  the  factor  for  stress  state  which  is 
most  important  in  the  development  of  internal  damage. 

The  initial  condition  for  isotropic  damage  is  uV(0)  =  0.  In  practical  analyses  the  coeffi¬ 
cient  v  is  of  the  order  10  (compare  ref.[4j).  Thus,  when  SI  units  are  used  in  the  analysis, 
the  factor  0  as  well  as  the  constant  H  reach  extremelv  high  values,  bevond  the  limit  of 
real  number  capacity  on  some  computers.  Thus,  for  numerical  analysis,  equation  (-1.12)  was 
recast  in  the  equivalent,  but  more  convenient  form: 


■J 


u 


where: 


Q  =  Q*  =  Aa+ax  -r  Bv/3X  +  Ci; 


H  =  (//-L 

sec 


The  additional  advantage  of  this  formulation  is  that  both  Q  and  H  are  in  the  stress  units 
(.UP,,)  instead  of  somewhat  cumbersome  ( MP^Y .  The  directional  damage  is  defined  in  a 
manner  very  similar  to  directional  hardening,  namely 


D  D 
:  =w  fjuu 


where  ii/j  are  directional  cosines  defined  in  equation  (4.11)  and  the  components  of  a  tensor 
u;0  evolve  according  to  equation 


D 

“W./ 


M 


where  </  and  M  are  material  constants.  The  initial  condition  is 


■fj( 0)  -  0 


.Vote  that  tliere  are  several  problems  with  practical  application  of  direct iona!  damage.  relia¬ 
bility  of  the  above  model  and  conducting  experiments  relevant  for  the  evolution  of  necessary 
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parameters.  Even  the  extensive  experiments  presented  in  references  [4.5.6]  did  not  provicU 
all  the  necessary  data  and.  hence,  the  damage  model  was  limited  to  the  isotropic  damage 
Therefore,  our  numerical  implementation  and  test  computations  were  also  limited  to  tin 
isotropic  damage. 


4.3  Equations  of  Equilibrium  and  Boundary  Conditions 


The  rate  version  of  the  equations  of  equilibrium  formulated  in  the  reference  configuration  is 


Div  (T)  +  b0  =  0 


(4.13 


where  T  is  the  first  Pioia  Kirchhoff  stress  tensor.  bQ  is  the  body  force  per  unit  volume  in  tin 
reference  configuration,  and  a  dot  denotes  time  derivative. 

After  substituting  the  relation  between  the  first  and  second  Piola-Kirchhoff  tensors  the 
above  equation  becomes 

Div  [FT)  +  b0  =  0  (4.14 


or  in  component  form 


j  +  b , 


=  0 


The  typical  boundary  conditions  are  prescribed  rates  of  displacements: 

ii  =  u{t)  in  r)f>mx[0.rj 


and  prescribed  rates  of  boundary  tractions 

TN=t0{t)  on  <9OH2x[0.r] 

or  in  terms  of  the  second  Piola-Kirchhoff  tensor 

(FT)N  =  t0(t)  on  09.R2  x  [0.  T\  (4.15 

In  the  above  Qr  is  the  domain  corresponding  to  the  reference  configuration  and  Ar  is  tin 
outward  normal  to  0ClR. 


4.4  Weak  Formulation 

In  order  to  obtain  a  weak  formulation  of  the  viscoplastic  evolution  problem  let  us  introduce 
the  space  of  test  functions: 


V  =  {v  £  [\V,n-p(nR)}\v  =  0  a.e.  on  rJQRl } 


and  the  space  of  trial  functions  (displacement  rates): 

L  =  (ii  G  [IKm-p(nn)]2,  u  =  u  a.e.  on  OP-m  } 

where  P.r  is  the  domain  occupied  by  the  reference  configuration  and  If  ’ /^ )  is  flic  Sobolev 
s])ace.  where  specific  values  of  m  and  p  depend  upon  the  particular  forms  of  constitutive 
ec|uations. 

Multiplying  the  equilibrium  equation  (4.14)  by  a  test  function  and  integrating  over  f !/? 
we  obtain  the  weak  form  of  the  rate  equilibrium  equations: 


Die  ( FT)  ■  vdX 


b  ■  vdx  Vtt  €  V 


The  application  of  the  divergence  theorem  to  the  left-hand  side  yields: 

[  (FT)  ■  I  V v )dX  =  [  60  •  vdX  +  [  (FT)N-vdS 
J  ri  pi  J  n  r  j  R 

and.  after  evaluation  of  the  time  derivative  on  the  left-hand  side,  we  obtain: 

/  (FT)  ■  (Vv)JX  ~  [  T  ■  (FTVv)dX  =  [  60  •  vdX  -f  [  (FT)N  ■  vdS 
JnR  J  nfi  JuR  JonR 

Now  it  is  easy  to  show,  using  basic  relations  presented  in  references  [15.23.36].  that: 

F  =  (/+Vu)  —  Vii 
and  by  symmetry  of  the  stress  tensor  T 

T  ■  ( Fr'Vv)  —  T  ■  6E(u.v) 

where  E  is  the  Green  strain  tensor  introduced  in  Section  3.1. 

Moreover,  we  have  from  equation  (4.15) 

(  FT)N  =  t0  on  r)QR2 

II  we  recall  the  decomposition  of  strain  rates  (4.S)  and  the  elastic  constitutive  equation- 
(4.9).  then  we  arrive  at  the  weak  formulation  of  the  evolution  problem: 

Find  a  displacement  rate  field  t — -u(x.t)  €  V  such  that 


/  V iiT  ■  V v  dX  +  /  CE  ■  6E(u.  v)  dX 

JnR  JnR 

=  f  CEin)  ■  6E(u.  v)  dX  +  /  CE(l)  ■  6E(u.  v ) 
J<.iR  J  n« 

+  /  b0  ■  v  dX  +  [  t  ■  v  dS  Vr  6  F 
J  Q  ft  J  dfl  ft 


(1.16) 
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With  the  index  notation,  this  equation  takes  the  form: 

/  «i .kTkjI\.J  dX  +  /  CuklEkL&Eij  dX 
JnR  JnR 


f  CuklE^SEu  dX  [  Cij,<lE{‘)l8E!j  dX 
JnR  JnR 


‘  i 


+  [  b0 ,  i’,  dX  +  f  toi^'i  db  Vi*,  £  V 
J Oh  J?Qr 

Note  that  for  the  large  displacement-small  strains  theory  we  have  det  F  =  1.  so  that  body 
forces  and  surface  tractions  can  be  referred  to  the  area  in  the  current  configuration: 

b0  —  b  in  f> 

t0  =  t  on  r)Q 


4.5  Finite  Element  Formulation 

A  finite  element  approximation  of  a  problem  (4.16)  is  obtained  by  introducing 

.v 

ulMX)  =  £u,a'Ia(X) 

a=i 

.V 

rf(X)  =^,A(X)  (4. IS) 

A  =  1 

«*(*)  =  £>,a'Pa(*) 

A  =  1 

where  .V  is  the  number  of  nodes.  'PA  is  a  shape  function  associated  with  node  A  and  :/,A 
is  the  value  of  uj1  at  node  A  (in  the  following  capital  Greek  symbols  will  be  reserved  for 
denumerating  nodes). 

I  lie  substitution  of  formulas  (4.18)  into  (4.17)  leads  to  the  following  matrix  equation 

(S  +  K)U  =  F{n)  +  F[t)  +  F{1)  (4. If)) 

Tl.  e  nonlinear  stiffness  matrix  in  each  element  can  be  calculated  as: 


'K=  GtCGJX 
J-q 


where  C  is  the  previously  introduced  elasticity  matrix  and  G  is  the  nonlinear  geometric 
matrix,  with  elem  ..its  defined  as 


C’UkA 


0E,j 

dllt cA 


dEu 

d"k  a 


IS 


The  occurence  of  the  matrix  S  in  the  above  equation  results  from  the  presence  of  the  defor¬ 
mation  gradient  F  in  the  rate  form  of  large  deformation  equilibrium  equation  (-1.14).  The 
elements  of  this  matrix  are  calculated  for  each  element  as: 

eSi±j  A  =  f  A.Kxb  a.lTkL  $1]  dX 

J'-Q 

and  represent  additional  stiffness  due  to  rotation  of  the  existing  current  stress  state. 

The  right-hand  side  vectors  represent  loads  due  to  the  rate  of  plastic  strains,  temperature 
and  mechanical  loads,  respectively.  These  vectors  are  defined  by 

eF(n]  =  [  GTCE[n)dX 

Tl,)  =  [GTElt)dX 

'F(,)  =  /  <PTb0dX  ~  [  $Tt0dS 
J'  O  Jd'-Q 

Note  that  it  is  convenient  to  represent  and  calculate  the  nonlinear  geometric  matrix  as: 


G  =  QQ 


or  in  component  form 

dFij  dum.x 

Gijk±  =  -X - T - 

OUm.X 

In  the  above.  C  is  a  purely  symbolic  geometric  matrix  obtained  by  differentiation  of  the 
nonlinear  kinetic  formula  (4.3).  while  Q  is  the  discretization  matrix  with  elements  defined 
by 

Q  m  A  k  ^  ^  k  m  ^  A .  X 

4.6  Viscoplastic  Solution  Method 

The  general  structure  of  the  finite  element  formulation  (4.19)  is  very  similar  to  that  obtained 
in  the  infinitesimal  theory  (see  Refs.  [33.35]).  the  difference  being  the  occurrence  of  the 
additional  matrix  5  and  the  dependence  of  matrix  K  on  the  current  deformation.  Thus  the 
general  computational  procedure  remains  practically  unchanged  after  the  extension  to  the 
large  displacement  theory. 

Also,  the  inclusion  of  a  damage  model  in  the  constitutive  equations  did  not  change 
the  general  structure  of  the  program,  but  did  require  additional  evolution  of  the  damage 
parameters. 

With  these  modifications,  the  new  algorithm  proceeds  through  the  following  steps: 
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1.  At  time  t.  initialize  T/j,  Zt,  and  aq  for  each  element. 

2.  Calculate  =  /^{Tu,  Z^,^' k)  for  each  element. 

3.  Assemble  and  solve  (5  -f-  K)U  =  F. 

4.  Calculate  Ejj  for  each  element.  E  —  GU. 

5.  Calculate  T{j  for  each  element  T  =  C(E  —  E (n))  —  CaT. 

6.  Calculate  Z,  for  each  element  Z,  =  g,(Tj Z,). 

7.  Calculate  A,  for  each  element  A,  =  h^Tix .  a;,-). 

S.  Integrate  u,-,  £j[)\  T/J-  Zt,  C,  forward  for  each  element  to  get  values  at  f  -r  At. 

9.  If  it  4-  AC  <  Cnai  g°  to  -•  otherwise  stop. 

In  the  case  of  the  predictor-corrector  method,  the  above  sequence  is  performed  at  both 
the  predictor  and  corrector  stages  and  the  final  values  are  obtained  using  averaged  rates. 

The  general  approach  for  this  scheme,  combined  with  the  adaptive  time  stepping,  is 
presented  in  Fig.  8.  Note  that  in  the  case  of  large  rotations  the  mechanical  deformation  can 
affect  the  solution  of  the  thermal  problem,  thus  providing  two-way  coupling  instead  of  the 
one-way  coupling  considered  so  far.  Analysis  of  such  cases  will  be  discussed  subsequently  in 
the  context  of  a  fully  coupled  fluid-thermal-structural  analyses. 

It  is  also  important  to  note  that,  due  to  the  application  of  a  rate  formulation,  extension  of 
the  analysis  to  large  displacements  did  not  cause  any  significant  increase  in  the  computational 
time:  in  particular,  no  additional  iteration  loops  were  necessary.  The  only  additional  cost 
was  a  slightly  more  complicated  evaluation  of  the  stiffness  matrix  and  load  vectors. 


5  Implicit/Explicit  Method  for  Compressible  Viscous 
Flow 

Finite  element  analysis  ot  compressible  flows  is  a  relatively  new  development,  but  siimificant 
progress  has  been  made  in  recent  years.  Computational  techniques  have  been  developed  for 
inviscid  and  viscous  flows.  Currently  there  exists  a  variety  of  successful  explicit  algorithms, 
the  most  popular  being  explicit  versions  of  the  Taylor- Galerh'in  schemes  [9. 12. 13. 27. 28. 29. 3lj 
and  Runge-Kutta  methods  [30.31].  Unfortunately,  as  indicated  by  our  numerical  experi¬ 
ments  and  the  experience  of  other  researchers,  explicit  algorithms  are  not  very  effective  in 
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the  analysis  of  viscous  Hows  and  in  the  prediction  of  aerothennal  loads  (pressure,  skin  fric¬ 
tion.  and  heat  flux).  In  these  problems,  due  to  stability  limitations  on  the  time  step,  they 
converge  very  slow.  As  a  result  current  research  in  this  area  is  focused  primarily  on  implicit 
schemes,  or  more  recently  implicit/explicit  schemes.  The  latter  approach  appears  to  be 
the  most  attractive  one.  as  it  combines  the  robust,  unconditionally  stable  implicit  schemes 
with  the  relatively  cheap  explicit  schemes,  to  achieve  maximum  effect iveness  at  minimum 
computational  cost. 

In  this  section  a  finite  element  formulation  of  an  adaptive  implicit/explicit  Taylor-Gal- 
erkin  method  for  compressible  viscous  flow  is  presented.  The  derivation  will  follow  several 
logical  steps.  First,  a  general  family  of  implicit  Taylor-Galerkin  algorithms  will  be  formu¬ 
lated.  Depending  on  the  selection  of  implicitness  parameters,  the  fully  explicit  scheme  or  a 
variety  of  implicit  schemes  can  be  obtained  within  this  family.  Secondly,  a  formulation  of 
an  implicit/ explicit  scheme  will  be  presented.  This  scheme  combines  implicit  and  explicit 
methods  in  different  subdomains  of  the  domain  Q.  Finally,  different  criteria  for  an  adap¬ 
tive  selection  of  implicit  and  explicit  zones  will  be  presented.  Several  numerical  tests  and 
validation  problems  will  be  discussed  in  Section  6.3. 

5.1  A  General  Family  of  Implicit  Taylor-Galerkin  Methods 

In  this  section  a  general  family  of  implicit  Taylor-Galerkin  methods  will  be  derived.  This 
family  is  based  on  a  combination  of  second-order  Taylor  series  expansions  in  time  with  a 
Galerkin  approximation  in  space  (one-,  two-,  or  three-dimensional).  Several  implicitness 
parameters  will  be  introduced,  so  that,  depending  on  the  particular  choice,  a  fully  explicit 
scheme  or  a  variety  of  implicit  schemes  can  be  recovered.  The  family  derived  here  includes, 
as  special  cases,  most  of  Taylor-Galerkin  type  algorithms  presented  in  the  literature. 

5.1.1  Basic  Derivation 
Navier-Stokes  Equations 

The  compressible  viscous  flow  of  a  calorically  perfect  gas  is  governed  by  the  Navier-Stokes 
erpiat ion  in  the  form: 

u  +  F,'t  =  F)'t  i  ■'>•!) 

where  u  is  the  vector  of  conservation  variables.  F,  are  i n viscid  (convective)  Dux  vectors, 
and  F\  are  viscous  flux  vectors.  Indices  i  in  the  above  formula  refer  to  axes  of  a  ('artesian 
coordinate  system,  a  comma  denotes  partial  differentiation  and  the  summation  convention 
is  applied.  The  components  of  vectors  u,  F,,  and  F ,l  are  given  in  two-dimensional 

u  =  {p.  pvi,  pv2,  pe}T  =  {p.  mi,  m>,  E}7 


raM*  hv: 


F,  = 


P'-'x 

PV\V,  +  p6u 

pVx  c\  4"  P&2 x 

( E  +  p)  v, 


F 


V 


0 

<7\x 

°~2, 

L'm  x  “t"  *7t 


(5.3) 

(5.-1) 


where  p  is  the  fluid  density,  p  is  the  pressure,  i\  are  velocity  components,  e  is  a  total  energy 
per  unit  mass,  <7,y  are  stress  components,  and  cp  are  components  of  a  heat  flux  vec  tor. 

It  should  be  noted  that  inviscid  fluxes  are  functions  of  the  conservation  variables  only: 


F ,  =  Ft(u) 


so  that  inviscid  Jacobians  are  defined  by: 


-4,  = 


dF\ 

du 


Viscous  fluxes  depend  on  both  conservation  variables  and  their  gradients: 


F]'  =  F)\u.u,} 


and  the  corresponding  Jacobians  are  defined  by 


Px 

RtJ 


f)Fx- 

du 


9F]' 

duj 


Second-Order  Taylor  Expansion  in  Time 

Assume  that  the  solution  un  is  given  at  the  time  moment  tn  and  the  solution  at  time  /'1  +  l 
is  to  l)e  calculated.  Formally,  the  values  of  the  solution  at  moments  tn  and  /'l_l  can  lie 
expressed  by  the  second -order  laylor  expansion  around  an  arbitral}'  moment  /  ~  '  imv  f  in. 
!)).  where  a  is  the  implicitness  parameter  with  values  between  zero  and  one: 

\/2 

un+1  =  un+d  +  (1  -  a)Atun+01  +  (1  -  a)2  — +  o(At3) 

At2 

u "  =  un+°  -  aA tun+°  +  a2  — un+c>  +  o(At 3) 


oo 


By  subtracting  these  two  formulas  one  obtains  a  formula  for  the  increment  oi  the  solution 
between  steps  n  and  n  +  1: 

A u  =  un+l  -  un  =  Atun+a  +  (1  -2a)  —  un+c>  +  of  At3)  (5.5^ 

Now  it  is  easy  to  observe  that: 

iin+°  =  un+0  +  o((q  -  3)At) 


so  that — still  preserving  second-order  accuracy — a  second  implicitness  parameter  can  be 
introduced  into  equation  (5.5): 

At2 

A u  =  Atun+a  +  (1  -  2a)  —  un+3  -f  o{AtJ)  (5.bi 

The  next  step  of  the  derivation  is  to  express  the  quantities  evaluated  at  time  moments  ’ 
and  tn *'1  by  quantities  evaluated  at  the  basic  steps  tn  and  tn+1.  It  is  easy  to  show,  using  a 
Taylor  series  expansion,  that: 

•un+l>  =  iin  4-  a  Ail  -f  o(At2) 
u'+3  =  iin  -T  3 Ail  +  ' i2) 

Substituting  these  formulas  into  equation  '"Tq  yields  a  two-parameter  expansion: 

At2 

Au  =  At(iin  4-  aAii)  4-  M  —  2o-)  —  (tzn  -f  3Au)  —  o(  At3)  <  5.7 ) 

Now.  following  the  original  idea  of  Lax  and  Wendroff.  the  original  equation  (5.1)  will  be  sub¬ 
stituted  to  ecfuation  (5.7)  to  replace  time  derivatives  by  space  derivatives.  This  substitution 
yields  a  formula  for  the  first  derivatives: 

u  =  F]'a  ~  F,.t  (5.5) 

and  for  the  second  order  derivatives,  after  substitution  of  the  inviscid  and  viscous  Jacobians: 


u  = 


a, (4 +  P, (fc  - F«) 


-  1.4.  F 


F, 


5.91 


It  can  be  observed  that  the  consecutive  terms  in  the  equation  (5.9)  involve  spatial  derivative's 
up  to  the  fourth  order.  Limiting  this  formula  to  terms  with  second-order  derivative's,  which 
can  be  effectively  handled  by  C°  continuous  finite  elements,  yields  an  approximation: 


u  =  [AiFk,k) A  +  o(/r,  k) 


(5.10) 
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where  o(fx,k)  represents  symbolically  a  quantity  of  the  order  of  viscosity  parameters  in  the 
Navier-Stokes  equations.  Substitution  of  formulas  (5.S)  and  (5.10)  into  the  incremental 
equations  (5.7)  gives  the  implicit  nonlinear  formula: 


An  =  At  [(JF*  -  *?.)  +  a  (*/£  -  AFU)] 

2  (5.11) 

+  (l-2a)—  [(AtJFfc.,),+ ^(A.F,,,),]  i-o(At3)+oiti.k)olAt2) 

It  can  be  noted  that  this  expansion  is  second-oruer  accurate  for  the  Euler  equations,  but 
only  first-order  accurate  for  viscous  problems.  In  the  latter  case,  for  the  sake  of  maximum 
generality,  a  third  implicitness  parameter  can  be  introduced  by  observing  that  AF ]  = 
k)o(At),  so  that  the  second  term  on  the  right-hand  side  of  equation  (5.11)  becomes 

aAF]t  —  qAF],  =  ~;AFUX  —  aAFx,x  +  (a  —  ~  )o(/.i.  k)o{At) 

This  substitution  preserves  the  first-order  accuracy  of  the  scheme. 

Linearization 

Equation  (5.7)  represents  a  nonlinear  formula  for  increments  of  the  solution  u  at  one  time 
step.  This  formula  is  nonlinear  due  to  nonlinear  dependence  of  the  fluxes  and  Jacobians  on 
the  solution  u.  It  can  be  shown,  however,  that  this  equation  can  be  effectively  linearized 
while  preserving  the  required  accuracy.  In  particular,  by  the  definition  of  the  inviscid  and 
viscous  Jacobians  one  can  write: 

AF,  =  A1-  Au  +  o(At2) 

AF]'  =  P"  Au,  +  R^Au.j  +o(At2) 

so  that  linearization  of  these  terms  preserves  second-order  accuracy  of  equation  (5.7b  Note 
that  if  the  dependence  of  viscous  fluxes  on  the  solution  u  is  neglected,  then 

AF]  =  RjjAUj  +  o(At) 


and  equation  (5.7)  is  only  first-order  accurate. 

The  second-order  term  in  the  equation  (5.7)  is  also  nonlinear  with  respect  to  u.  due  to 
dependence  of  A,  on  the  solution.  Since,  however,  this  term  is  multiplied  by  At2,  it  actually 
suffices  to  observe  that  =  o{At)  and  use  an  approximate  expansion: 

M  A,Fk.k).t  =  ( A?AnkAu  fc),  +  o(At) 
in  order  to  preserve  the  second-order  accuracy  of  the  basic  formula  (5.7). 


After  substitution  of  these  linearized  terms  into  the  equation  (5.7)  and  regrouping,  a 
general  linear  implicit  incremental  formula  is  obtained: 


Au  +  a  At  ( A"  Au)  ^  —  7  At  (R^Au.j')  .  +  (P^Au)t 


+ 


-  (1  -2a)^(A;A^uj). 


=  A (  (Fp  -  FI,)  +  (1  -  2a)44  (AfFlj, 


(5.12) 


The  above  equation  is  combined  with  the  appropriate  boundary  conditions  to  form  a  bound¬ 
ary-value  problem  over  the  domain  fi.  The  particular  forms  of  boundary  conditions  and 
their  finite  element  implementations  are  discussed  in  Section  5.4.  The  incremental  solutions 
at  consecutive  time  steps  are  accumulated  to  formulate  a  solution  of  the  time-dependent 
problem.  If  one  prefers  to  do  so.  the  above  equation  can  easily  be  recast  in  a  form  in  which 
the  new  solution  un+1  is  calculated  instead  of  the  increment  Au. 


Remarks  Concerning  Eq.  (5.12) 

1.  Formula  (5.12)  is  a  linear  equation  for  the  increment  Au.  The  formulation  is  in  general 
implicit,  except  when  a  =  , 3  =  7  =  0,  in  which  case  the  explicit  scheme  is  recovered. 
The  interpretation  of  different  implicitness  parameters  is  as  follows: 

a  —  represents  implicitness  of  the  advective  terms. 

3  —  represents  implicitness  of  the  viscous  terms. 

7  —  represents  implicitness  of  the  second  order  terms. 

2.  The  scheme  is  generally  first-order  accurate  in  time  for  the  Xavier-Stokes  equations, 
but  it  is  second-order  accurate  for: 

—  Euler  equations 

—  Xavier-Stokes  equations  with  0  =  7  =  .5  (Crank-Xicholson  scheme) 

In  order  to  distinguish  scheme  (5.12)  from  the  purely  first-order  methods,  it  will  be 
referred  to  as  an  (incomplete)  second-order  scheme. 

•7  If  all  the  terms  with  coefficient  (1  —  2a)  are  eliminated  in  equation  (5.12).  then  the 
purely  first-order  scheme  is  obtained. 

1.  The  above  scheme  includes  as  special  cases  most  of  the  other  finite  element  schemes 
based  on  a  second-order  Taylor  series  expansion  in  time,  in  particular: 
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•  a  =  0,  ,3  =  0.  7  =  0  corresponds  to  the  popular  explicit  Taylor-Galerkin  method. 

•  a  =  0.  /3  ^  0  with  no  viscosity  terms  corresponds  to  the  so-called  weakly  implicit 
method  applied  to  the  Euler  equations  by  Demkowicz.  Rachowicz.  and  Oden  [11] 
and  others, 

•  q  =  0,  $  =  1,  7  =  1  corresponds  to  an  implicit  scheme  used  by  Hassan.  Morgan, 
and  Peraire  [17]. 

Moreover,  a  variety  of  second-order  algorithms  developed  by  Donea  [9]  and  Donea. 
Giuliani,  Laval,  and  Quartapelle  [12]  can  be  obtained  within  this  family  (with  the 
exception  of  the  third-order  schemes  developed  by  Donea.  Quartapelle  and  Selmin  [13] 
for  convection  problems). 

Similarly,  several  corresponding  finite  difference  algorithms  can  also  be  identified.  In  partic¬ 
ular  a  family  of  implicit  schemes  developed  for  the  Euler  equations  by  Lerat  [21]. 

5.1.2  Linear  Stability  Analysis 

It  is  of  interest  to  analyze  stability,  dissipation,  and  dispersion  properties  of  the  various 
schemes  obtained  within  the  family  derived  in  the  previous  section.  Such  an  analysis  was 
performed  for  a  simplified,  one-dimensional  linear  advection-diffusion  equation  of  the  form: 

u  -f  au  x  +  b\i'XX  —  0  in  [0,  L]  x  [0.  T] 

with  periodic  boundary  conditions: 

u(0.t)  =  u(L.t) 

MO.*)  =  u.x{L.  t) 

For  this  equation,  incremental  equation  corresponding  to  the  formula  (5.121  was  formulated 
and  a  discrete  finite  element  scheme  was  derived,  assuming: 

-  a  uniform  mesh  of  spacing  h  —  L/ . V.  where  X  is  the  total  number  of  elements 

-  a  lumped  mass  matrix 

The  lumped  mass  matrix  was  considered  here  because — as  results  from  arguments  which  will 
be  presented  in  Section  5.4 — it  significantly  improves  the  computational  efficiency  of  mixed 
implicit/explicit  schemes.  Application  of  a  standard  von  Neumann  stability  analysis  to  the 
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resulting  scheme,  yields  a  formula  for  an  amplification  coefficient  for  mode  m(m  =  0 
in  the  form: 


-V) 


c  =  - 

^  m  — 


1  + 

C  F  L  1 

((1  -7)~5— +  (1  -2a)(l  -  3)C F L2)(cos  3m  -  1) 
lie 

—  u(  1  —  a  1C F L  sin  T,„] 

1  +  (  7  Re  +  (*  3)CFLSj  (cos  3m  1) 

-r  t[—aCFL  sin  3,„  ] 

(o. 


where  CFL  is  a  Courant-Friedrich-Levy  number  for  this  problem: 

aAt 


Re  is  a  mesh  Revnolds  number: 


CFL  = 


ah 

Re=T 


3m  is  defined  in  terms  of  a  wave  number  km  as: 


3m  —  km  h  — 


mh 

~L 


and  i  =  \/— T- 

The  above  formula  was  used  to  analyze  the  dissipation  and  dispersion  properties  for 
various  combinations  of  implicitness  parameters.  Some  results  of  the  stability  analysis  are 
summarized  below.  For  simplicity,  the  figures  presented  will  refer  only  to  pure  advection 
problems  (1  /Re  =  0). 


1.  Fully  explicit  scheme:  a  =  0.  3  =  0.  *  =  0. 

The  results  of  stability  analysis  for  this  scheme  are  graphically  presented  in  Fig.  10. 
The  figure  presents  modules  of  the  amplication  factor  C  for  increasing  values  of  the 
CFL  number  (which  corresponds  to  an  increasing  time  step).  Actually,  the  figure  is 
a  graph,  with  values  of  the  amplification  factors  for  different  modes  spanning  the  area 
between  a  horizontal  line  |Cj  =  1  (mode  m  —  0).  and  curve  4  (mode  m  —  A  ).  The 
scheme  is  stable  for  given  time  step  if  the  graph  for  this  step  is  below  the  horizontal 
line  |C'|  =  1.  otherwise  the  scheme  is  unstable. 

As  expected,  the  explicit  scheme  is  conditionally  stable  with  the  stability  limit  for  pure 
advection  being  CFL  <  1. 

2.  Weakly  implicit  scheme:  a  =  0,  3  =  .5.  7  =  .5. 

The  stability  and  dissipation  analysis  for  this  scheme  is  presented  in  Fig.  11.  The 
scheme  is  unconditionally  stable,  with  dissipation  properties  diminishing  with  an  in¬ 
creasing  time  step.  Thus  it  can  be  expected  that  this  scheme  may  be  oscillatory  at 
large  time  steps. 


27 


3.  Weakly  implicit  scheme:  a  =  0,  ,5  =  1,  7  =  1. 

The  stability  and  dissipation  analysis  for  this  scheme  is  presented  in  Fig.  12.  The 
scheme  is  unconditionally  stable,  with  dissipation  properties  for  higher  modes  increas¬ 
ing  with  an  increasing  time  step. 

4.  Crank-Nicholson  scheme:  a  =  .5,  $  =  .5.  7  =  .5. 

This  scheme  is  unconditionally  stable  (See  Fig.  13),  with  no  dissipation  for  pure 
advection  problems  (for  a  regular  one-dimensional  mesh).  Thus  it  can  be  expected 
to  be  oscillatory  for  advection-dominated  problems.  However,  since  this  scheme  is 
actually  second-order  accurate  for  the  Navier-Stokes  equations,  it  can  be  very  effective 
within  boundary  layer  subregions. 

Before  other  schemes  are  discussed,  it  should  be  noted  that  for  the  schemes  discussed  so  far. 
the  implicitness  parameter  a  was  equal  to  zero.  The  following  two  cases  will  have  o  A  0. 

•5.  Implicit  scheme:  a  =  1.  3  =  0.  7  =  0. 

The  stability  and  dissipation  analysis  for  this  scheme  is  presented  in  Fig.  14.  Surpris¬ 
ingly,  implicit  treatment  of  the  convection  terms  spoils  the  stability  and  this  scheme 
is  unconditionally  unstable  for  pure  advection  problems. 

6.  Fully  implicit  scheme:  a- =  1.  ,5=1,  7  =  1. 

Contradictory  to  what  would  be  intuitively  expected,  this  fully  implicit  scheme  is  not 
the  most  stable.  As  can  be  seen  in  Fig.  1-5.  this  scheme  is  stable  for  large  time  steps 
and  unstable  for  small  time  steps.  This  surprising  result  is  actually  confirmed  in  a 
numerical  solution  of  the  Navier-Stokes  equations. 

The  above  few  examples  outline  the  nature  of  dependence  of  the  stability  and  dissi¬ 
pation  properties  of  second-order  schemes  on  the  selection  of  the  implicitness  parameters. 
The  behavior  of  schemes  corresponding  to  intermediate  values  of  these  parameters  can  be 
interpolated  from  the  above  results  or  precisely  analyzed  using  the  formula  (o.l3). 

A  similar  stability  analysis  can  be  performed  for  the  first-order  scheme,  which  is  obtained 
by  dropping  terms  with  multiplier  (1  -2a).  In  this  case,  the  coefficient  .5  disappears  from 
equations  and  properties  of  the  scheme  depend  on  the  values  of  a  and  7.  The  general  result 
of  the  stability  analysis  is  that: 

•  For  a  <  .5  and  7  <  .5  the  scheme  is  conditionally  stable  for  advection-diffusion 
problems.  For  pure  advection  problems,  t he  scheme  is  unconditionally  unstable. 
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•  For  .5  <  a  <  1  and  5  <  7  <  1  the  scheme  is  unconditionally  stable.  The  dissipation 
and  dispersion  properties  for  a  =  1,  7  =  1  are  quite  similar  to  the  properties  of  the 
second  order  scheme  with  a  =  0.  3  =  1.  7  =  1. 

It  is  of  interest  to  note  that  the  effect  of  an  implicit  treatment  of  the  convective  terms 
changes  dramatically  depending  on  the  order  of  the  scheme.  While  in  the  first-order  scheme 
it  is  necessary  for  unconditional  stability,  in  the  second-order  scheme  it  actually  destroys 
stability. 

The  results  of  the  above  stability  analysis,  although  based  on  the  simplified,  one-dimen¬ 
sional  equation,  are  very  well  confirmed  in  the  solution  of  the  Navier- Stokes  equations  in 
multidimensional  domains. 

5.2  Weak  Formulation 

In  order  to  obtain  a  weak,  variational  formulation  of  the  incremental  equation  (5.12).  the 
spaces  of  trial  and  test  functions  were  introduced: 

U  =  {u  =  ( t/ ] ,  u2 . . .  u\i)  s.t.  u,  €  Hl{Q)  and  u,  =  u,  on  I'd,} 

I  =  {t<  =  (ipri.-.ru)  s.t.  q  6  and  r,  =  0  on  F/j,} 

where  M  is  the  number  of  conservation  variables.  Hl{Cl)  is  the  Hilbert  space,  and  I'o,  is 
the  boundary  with  specified  Direchlet  boundary  conditions  for  variable  it,.  Note  that  the 
notion  of  Dirichlet  boundary  conditions  is  rather  symbolic  here,  since  the  actual  boundary 
conditions  are  often  imposed  on  characteristic  variables  rather  than  conservation  variables 
a,  (see  Section  5.41. 

After  multiplication  of  the  incremental  equation  (5.11)  by  an  arbitrary  test  function  iqxl. 
integrating  over  the  domain  H  and  application  of  the  divergence  theorem,  the  following  weak 
formulation  of  the  problem  is  obtained: 
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=  ^  A t  (f:  -  F'f")  ■  ® -  (1  -  2 ^  —  A^F^  •  VjdSl 

4-  J  —Ain,  (>"  -  *Tn)  •  »  +  (1  -  2a)=^-niA?FlJ  ■  vdS 

It  can  be  shown,  by  selection  of  appropriate  test  functions,  that  the  solution  of  this  prob¬ 
lem  is — in  the  sense  of  distributions — the  solution  of  the  boundary-value  problem  formulated 
in  Section  5.1. 


5.3  Finite  Element  Approximation 

The  Galerkin  approximation  of  the  variational  problem  (5.14)  is  obtained  by  partitioning  the 
domain  Q  into  the  finite  element  mesh  and  introducing  an  approximation  of  trial  functions 
and  their  derivatives: 

.v 

uh(x.t)  =  YlU,*i(x) 
i~  i 
.V 

u^x.t)  =  TTUi'VkU r) 

/=  i 

and  the  same  for  test  functions  v.  In  the  above  formula  .V  is  the  number  of  nodes  within 
the  domain  Q  and  vI*/(x)  is  the  shape  function  associated  with  node  I.  After  introducing 
this  approximation  into  the  variational  statement  (5.14)  and  requiring  that  it  be  satislieu 
lor  every  combination  oi  V/.  the  following  matrix  equation  is  obtained: 

(M  +  K)AU  =  R 

where  M  is  the  mass  matrix,  K  is  the  ''stiffness”  matrix  and  R  is  the  right-hand  side  vector. 
It  is  convenient  to  present  the  above  equation  in  the  semi-indexed  notation: 

{Mu  +  K  u)AU  j  =  R[ 
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where  Mu  and  K\j  are  small  blocks  of  size  M  x  M  associated  with  nodes  /  and  ./ 
that  M  is  the  number  of  conservation  variables). 

The  particular  forms  of  these  matrices  corresponding  to  the  variational  problem 
is — before  application  of  boundary  conditions — the  following: 

1.  Mass  matrix 

Mu  =  f  /.ux.u'I'  1$ jdP- 

Jn 

2.  Stiffness  matrix  due  to  interior  integrals: 

K,j  =  [  -aAtA^j'Vu 

Jq 


+ 


+  7  A 

+  (1  —  2a)3^2-A™  A*xVj'jxV[',dQ 

3.  Right-hand  side  due  to  interior  integrals: 

R,  =  jn  M  (F:  -  Ftr"  )  . 

-  (i 

4.  Stiffness  matrix  due  to  boundary  integrals: 

Ku  -  f  ou±tA,nt'Vj'$[  + 

Jen 


-  •■AtPlnt'Vj'Vi  4- 

At2 

-  ( 1  -  2a) p  —  n,A,Aj  T  j. :  T  / ds 

5.  Right-hand  side  due  to  boundary  integrals: 

R,  =  [  -AtFlmVi 

Jdn 


( recall 
(5.14) 
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+  A  tF)'nniVr 


+  {l-2a)—niA?F]'J*t 

It  should  be  noted  that — as  a  result  of  the  derivation  of  the  algorithm — the  elements  of 
the  stiffness  matrix  K  are  in  fact  consistent  linearizations  of  corresponding  terms  on  the 
right-hand  side,  for  example: 

Kv  .JR] 

where  superscript  V’  denotes  terms  associated  with  natural  viscosity.  This  remark  is  very 
important  for  the  development  of  various  modifications  of  the  algorithm,  because  every 
modification  of  the  right-hand  side  must  be  accompanied  by  a  consistent  modification  to 
the  left-hand  side — otherwise  the  stability  of  the  algorithm  is  spoiled.  For  example,  if  the 
group  approximation  of  convective  fluxes  is  applied,  then  the  calculation  of  .Jacobians  .4,  on 
the  left-hand  side  must  correspond  to  this  procedure. 

5.4  Artificial  Dissipation 

In  order  to  suppress  spurious  oscillations  of  the  solution,  a  variety  of  models  of  artificial 
dissipation  are  used.  In  this  work  we  will  assume  that  the  artificial  dissipation  can  be 
introduced  as  the  additional  flux  in  the  Xavier-Stokes  equations  in  the  form: 

u  +  F„  =  F]  t  +  F;4, 

where  the  artificial  dissipation  flux  is  the  function  of  the  solution  vector  and  its  derivatives: 

F;4  =  F;V.u,) 


with  corresponding  Jacobians  defined  as: 


P;4  = 


The  advantage  of  this  approach  is  that  the  artificial  dissipation  can  be  treated  using 
exactly  the  same  formulation  and  procedures  as  for  the  natural  viscosity.  In  the  implicit 
algorithm,  for  the  sake  of  generality,  a  fourth  implicitness  parameter  7  is  introduced  for  the 
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terms  associated  with  the  artificial  dissipation.  In  the  calculation  of  the  stiffness  matrices, 
right-hand  sides  and  boundary  terms  the  same  formulas  are  used  as  for  the  natural  viscosity. 

Within  the  above  framework,  various  models  of  artificial  dissipation  can  be  formulated 
relatively  easy.  For  example,  a  straightforward  extension  of  the  original  Lapidus  dissipation 
[20]  to  multidimensional  cases  yields  artificial  fluxes  defined  as: 


F;4  =  F'uj 


(5.15) 


k  —  CkAe  [U|t|  j 

where  c ^  is  a  coefficient  (usually  between  zero  and  one),  Ae  is  the  element  area,  and  c,  are 
the  components  of  velocity  vector  (no  summation  on  i).  The  Jacobians  P  and  R  can  be 
defined  by  straightforward  differentiation  of  formula  (5.15). 

Another  generalization  of  the  Lapidus  concept  was  proposed  by  Lohner  and  Morgan  [22]. 
Within  the  framework  proposed  here,  the  fluxes  corresponding  to  their  model  are  of  the 
form: 

or 

F;4  =  likljUj  (5. 16) 

where  /  is  the  normalized  gradient  of  the  magnitude  of  velocity: 

l  =  9rad\v\ 

|(jrrad|v|| 

and  the  coefficient  k  is  calculated  from  the  formula: 


k  =  ciAAl  ■  nrad{v  •  /)) 

The  .Jacobians  P  and  R  can  be  defined  by  differentiation  of  formula  (5.15).  If.  for  simplicity, 
dependence  of  k  and  l  on  the  solution  is  disregarded,  then: 

Pi  =  0 

(5.17) 

R,j  =  kljjl 

where  I  is  the  identity  matrix  of  dimension  U.  In  the  incremental  equation  (5.12).  the  above 
approach  leads  to  an  additional  term  of  the  form: 

,  ,  ,  k)u 

l,k— 

al 
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which  differs  slightly  from  the  original  formula  proposed  by  Lohner  and  Morgan  (equation 
(9)  in  [22]).  These  two  versions  are  equivalent  only  for  /  constant  throughout  the  domain. 
It  can  be  verified,  however,  that  the  directional  derivative  used  in  [22]  is  not  in  divergence 
form  and  thus  cannot  be  directly  used  in  the  variational  formulation  for  arbitrary  /.  It  is 
actually  the  above  term  that  yields — with  no  additional  assumptions — the  diffusion  matrix 
R  of  the  form  (5.17). 

5.5  Boundary  Conditions 

Proper  formulation  and  application  of  boundary  conditions  is  one  of  the  most  difficult  prob¬ 
lems  of  computational  fluid  mechanics.  While  this  topic  has  received  considerable  attention 
from  the  theoretical  point  of  view  and  within  the  finite  difference  method  framework  (ref¬ 
erences  [14.16.32]).  a  theoretical  background  for  the  treatment  of  the  boundary  conditions 
within  the  finite  element  framework,  especially  in  an  implicit  approach,  is  not  yet  complete. 

Before  discussion  of  boundary  conditions  for  the  implicit  Taylor-Galerkin  method,  it  is 
useful  to  quote  the  general  result  of  Strikverda  [32],  which  specifies  the  number  of  boundary 
conditions  necessary  for  well-posedness  of  the  linearized  Euler  and  Navier-Stokes  equations. 
These  results  are  summarized  for  two-dimensional  problems  in  the  table  below. 


Type  of 

Boundary 

Euler 

(not  regularized) 

Xavier-Stokes 

supersonic  inflow 

4  ess 

4  ess 

subsonic  inflow 

3  ess 

3  ess  4-  1  nat 

subsonic  outflow 

1  ess 

1  ess  -f  2  nat 

supersonic  outflow 

0 

0  ess  +  3  nat 

no-flow 

I  ess 

1  ess  +  2  nat 

solid  wall 

— isothermal 

— 

3  ess 

— heat  flux 

_ 

2  ess  +  1  nat 

In  this  table  "ess"  denotes  the  essential  boundary  condition  and  "nat"  denotes  natu¬ 
ral  boundary  conditions.  The  essential  conditions  are  to  be  imposed  on  the  characterist  ic 
variables  rather  than  the  conservation  variables  (see  Section  5.5.1).  It  is  of  importance  to 
note  that  the  numbers  presented  in  the  table  are  true  for  problems  that  are  not  regularized. 
If — as  is  the  case  of  virtually  all  computational  techniques — some  artificial  diffusion  is  built 
into  the  algorithm  or  added  explicitly,  natural  boundary  conditions  should  be  imposed  on 
these  terms  even  for  Euler  problems.  Moreover,  since  artificial  diffusion — in  contradiction  to 
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the  natural  viscosity — affects  all  the  conservation  variables,  the  number  of  natural  bound¬ 
ary  conditions  for  these  terms  should  actually  be  one  more  than  for  the  (nonregularized) 
Xavier-Stokes  equations. 

In  the  following  sections  the  formulation  of  the  boundary  conditions  for  various  boundary 
types  will  be  discussed.  Presentation  of  the  theoretical  background  will  be  limited  to  the 
necessary  minimum  and  the  interested  reader  is  referred  to  works  [14.16.32]  and  references 
therein. 

5.5.1  Characteristic  Decomposition  of  Boundary  Terms 

Before  the  presentation  of  the  various  boundary  conditions  it  is  useful  to  recast  the  boundary 
integrals  in  terms  of  the  characteristic  variables  (Riemann  invariants). 

In  the  variational  formulation  a  typical  boundary  term  is  of  the  form: 

AxntAu  ■  v  =  AnAu  ■  v  (5.18) 

where  .4n  is  a  nonsymmetric  matrix.  The  matrix  An  can  formally  be  presented  in  its  own 
eigenbasis  as: 

M 

An  —  'y  '  ^a{bo,  S  C.y) 

<1=  1 

where  6L>  and  c.y  are  the  left  and  right  eigenvectors,  respectively.  The  eigenvalues  A.,  are 
defined  in  the  two-dimensional  case  as: 

Ai  =  vn  —  c 

A?  =  l'n 

A3  =  Vn 

A  4  =  Vn  +  C 

where  v„  is  the  velocity  normal  to  the  boundary  and  c  is  the  speed  of  sound.  The  expressions 
for  eigenvectors  6a  and  c,y  can  be  found  in  references  [11.14.32].  Note  that  a  positive  value 
of  \0  means  that  the  corresponding  characteristic  exits  the  domain  across  a  given  boundary, 
while  negative  values  of  A  ,  correspond  to  signals  entering  the  domain. 

I  lie  characteristic:  variable  ig  are  defined  as  components  of  the1  vector  u  in  the  eigenbasis 
of  .4n  so  that: 

Au  =  (Au  ■  6a)ca  =  A uac0 
V  =  (v  ■  ca)ba  =  vab0 
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With  these  definitions,  the  boundary  formula  (5. IS)  can  be  presented  in  terms  of  character¬ 
istic  variables  as: 

M 

AnAu  ■  v  =  (5.1!)) 

i>  =  l 

The  above  representation  is  very  useful  in  the  formulation  of  essential  boundary  conditions. 
As  a  general  rule,  the  characteristic  variables  corresponding  to  characteristics  entering  the 
domain  (negative  Aa)  need  to  be  specified  as  the  essential  boundary  conditions,  while  the 
characteristic  variables  exiting  the  domain  are  continued  across  the  boundary  from  the  in¬ 
terior. 


5.5.2  Supersonic  Inflow 

On  the  supersonic  inflow  the  values  of  all  the  characteristic  variables  (thus  also  of  all  the 
conservation  variables)  are  specified  as  the  upstream  values.  Formally,  this  means  that: 

u  =  it  on  OCl  i 


or.  in  the  incremental  form: 

iu  =  u  —  un  on  dO. i 

In  practical  applications  these  conditions  are  enforced  by  the  penalty  method,  which  is 
obtained  by  adding  to  the  variational  formulation  a  term: 


1 

€ 


—  (u  ~  un )]  •  vds 


where  e  is  a  small  parameter. 

The  corresponding  stiffness  matrices  and  right-hand  sides  are  of  the  form: 


Ku  =  /  -Nttyjds 

Jan,  e 

Ri  =  /  -(u-unVVjds 

Jan,  e 

where  I  is  the  identity  matrix  of  the  size  M .  Simple  nodewise  enforcement  of  these  conditions 
is  obtained  by  formally  replacing  the  shape  functions  with  Dirac  delta  functions  t<>  obtain: 


I<u  =  -f 

e 

R,  =  ~{U i  -  Un{) 

t 
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5.5.3  Supersonic  Outflow 


On  the  supersonic  outflow  all  the  characteristic  variables  are  continued  from  the  interior  of 
the  domain,  so  that  no  essential  boundary  conditions  are  imposed.  The  explicit  algorithm 
application  of  supersonic  outflow  boundary  conditions  is  achieved  simply  by  calculation  of 
the  boundary  integrals.  In  the  implicit  algorithm,  natural  boundary  conditions  are  imposed 
on  viscous  and  second-order  terms. 

Natural  boundary  conditions  on  viscous  terms  are  imposed  by  observing  that  the  viscous 
boundary  terms  on  the  left-hand  side  of  the  variational  equation  (5.14)  can  be  interpreted 
as: 


K]jAUj  =  -7A<(JR"nt^'I'/  +  P1n,'P;VlrJ)  AUjds 

=  /  —.AtAF^Vids 

J  dQ0 

where  the  components  of  AF\  are: 

AF\  =  {0.  Acrln,N(J,n,  Aqn}7 

In  order  to  formally  impose  natural  boundary  conditions  the  above  terms  are  transferred  to 
the  right-hand  side  with  prescribed  values  of  AF\  .  so  that  the  new  right-hand  side  is: 

R[  =  R[  +  /  7 AtAF  T  [ds 

Jen* 

Note  that  since  the  mass  flux  due  to  viscous  terms  is  identically  zero,  this  procedure  actually 
imposes  only  three  natural  boundary  conditions  (in  two  dimensions). 

ine  choice  of  the  actual  conditions  is  somewhat  arbitrary.  Currently  two  options  are 
implemented,  namely: 

•  zero  change  of  flux  at  the  time  step  (frozen  viscous  flux): 

AFl=o 

•  zero  total  viscous  llux.  enforced  by: 

3f:: = o  - 

Obviously,  on  an  outflow  boundary  adjacent  to  a  solid  wall  the  viscous  fluxes  arc'  not  zero 
and  the  first  procedure  is  more  appropriate. 
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In  order  to  ensure  well-posedness  of  the  problem,  proper  natural  boundary  conditions 
should  also  be  imposed  on  the  second-order  terms.  Analogously,  as  in  the  viscous  case, 
second-order  terms  on  the  left-hand  side  of  equation  (5.14)  can  be  transformed  to  the  form: 

KfjAUj  =  [  -(l-2a)0—-^I(AnAF,.J)ds 

van,  2 

This  term  has  no  simple  physical  interpretation  and  therefore  the  selection  of  the  natu¬ 
ral  boundary  condition  to  be  imposed  is  somewhat  difficult.  The  procedure  adopted  by 
Demkowicz.Rachowicz  and  Oden  [11]  is  to  decompose  the  above  term  into  components  nor¬ 
mal  and  tangential  to  the  boundary  and  impose  boundary  conditions  only  on  the  normal 
term  (symmetry  boundary  conditions).  In  this  work  a  slightly  different  procedure  was  ap¬ 
plied.  according  to  which  the  whole  term  is  transferred  to  the  right-hand  side  with  certain 
prescribed  values.  This  corresponds  to  imposing  natural  boundary  conditions  on 

AnAFhJ  =  AnAjAuj 

The  actual  boundary  condition  applied  is  to  set  total  value  of  this  term  to  zero,  so  that  the 
prescribed  value  at  the  time  step  is: 

AnAjAuj  =  0  -  AnAjUn:  (5.20) 

The  intuitive  interpretation  of  this  condition  can  be  obtained  by  observing  that  for  Euler 
problems  the  enforced  condition  is  Anii  =  0.  or  in  terms  of  characteristic  variables: 

\ailabc  =  0  Q  =  1 . M 

Since  on  the  supersonic  outflow  all  the  eigenvalues  satisfy  A*  >  0.  this  condition  means 
that  the  characteristic  variables  do  not  change  in  time  as  they  cross  the  boundary.  This 
corresponds  quite  well  to  the  idea  of  the  continuation  from  the  interior  across  the  far— field 
boundary. 

A  somewhat  more  appealing  interpretation  can  be  presented  for  the  simple  two-dimen¬ 
sional  advection  equation: 

;’/  =  (1,11,  i  =  1.2 

for  which  tlu*  character ist  ics  are  straight  lines  defined  in  space-time  bv  the  vector  c  = 

{ a ! ,  cr 2,  1 }  (see  Fig.  Lb).  The  natural  boundary  condition  corresponding  to  (5.20)  is: 

nta,ajU  :  =  0 


or  equivalently: 


an(Du.  a)  =  0 


3S 


where  (Du.  a)  denotes  directional  derivative  of  u  in  the  direction  of  a.  This  means  that  on 
the  outflow  boundary  (an  >  0)  u  must  be  constant  along  advection  lines. 

On  the  contour  plot  this  enforces  contours  of  u  to  be  parallel  to  advection  lines  on  the 
boundary.  It  can  be  observed  that  the  condition  applied  by  Demkowicz.  Rachowicz  and 
Oden  [II]  enforces  derivatives  of  u  along  the  normal  to  the  boundary  n  to  be  zero,  which 
corresponds  to  contour  lines  normal  to  the  boundary. 

5.5.4  Subsonic  Inflow  and  Outflow 

On  the  subsonic  inflow  or  outflow  boundary  essential  boundary  conditions  are  imposed  only 
on  the  characteristic  variables  corresponding  to  negative  eigenvalues  A.-,.  For  each  of  these 
terms  the  corresponding  condition  is: 

ua  =  u.y  on  dO.io 


or,  in  the  incremental  form: 

A ua  =  uQ  -  u* 

If  the  prescribed  far-field  values  of  the  conservation  variables  are  denoted  by  u.  then  the 
penalty  term  enforcing  essential  boundary  conditions  on  the  selected  characteristic  variables 
uQ  is  of  the  form: 

Ku  =  -/  (c*  ®  ba)'V  1$ jds 

Rj  =  -  /  [(u-un)-6a]<VM$ 

£  J  dU 

Xodewise  application  of  these  conditions  is  obtained  analogously  as  for  supersonic  inflow  by 
replacing  the  shape  functions  with  Dirac  delta  functions. 

For  the  characteristic  variables  with  nonnegative  eigenvalues  the  "continuation  from  the 
inside”  conditions  should  be  applied.  Formulating  these  conditions  for  selected  characteristic 
variables  leads  to  rather  complicated  formulas.  Thus,  for  practical  applications  it  is  better 
to  observe  that,  since  the  penalty  procedure  actually  overrides  any  other  conditions,  the 
continuation  condition  can  be  first  applied  to  the  whole  vector  of  conservation  variables  (by 
the  supersonic  outflow  procedure),  and  then  the  above  penalty  method  can  be'  applied  to 
selectively  enforce  essential  boundary  conditions. 

It  should  also  be  noted  that  all  the  farfield  boundaries  are  actually  artificial  and  are 
supposed  to  represent  the  cut-oil  sections  of  the  flowfield  where  "nothing  special  happens." 
Thus,  if  these  boundaries  are  chosen  too  close  to  the  complex  flow  phenomena,  ambiguous 
results  can  always  be  expected. 
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5.5.5  Solid  Wall 


There  are  basically  two  types  of  solid  wall  boundaries: 

•  adiabatic  walls  with  prescribed  zero  heat  flux  (.YZ-th  component  of  the  viscous  flux 
vector): 

9n  =  Fn(M)  ~  0 

•  isothermal  walls  with  specified  temperature: 


T  =  T 


In  addition  to  the  above  conditions,  zero  velocity  (zero  momentum)  conditions  are  also 
specified  on  the  solid  wall.  These  conditions  are  easily  enforced  by  the  selective  application 
of  the  penalty  method,  similar  to  the  supersonic  inflow  procedure.  The  incremental  form  of 
the  adiabatic  condition  of  zero  heat  flux  is 


AF 


vn 

n(A/) 


(5.21) 


This  natural  boundary  condition  is  applied  by  formally  transferring  the  viscous  terms  cor¬ 
responding  to  the  energy  equation  from  the  left-hand  side  of  the  variational  equation  (5.1-1) 
to  the  right-hand  side  and  setting  the  increment  of  the  heat  flux  according  to  (5.21).  This 
type  of  procedure  was  already  presented  for  supersonic  outflow.  It  is  of  interest  to  note  that 
since  the  viscous  terms  do  not  directly  affect  mass  fluxes  and  the  momentum  equations  are 
overridden  by  the  penalty  method,  in  practice  all  viscous  contributions  can  be  skipped  on 
the  left-hand  side. 

On  the  isothermal  wall  the  additional  boundary  condition  is  a  prescribed  temperature  T 
or.  equivalently,  a  prescribed  specific  energy  e.  Since  the  kinetic  energy  is  zero  on  the  wall, 
the  above  condition  can  be  expressed  in  terms  of  conservation  variables  as: 

E  _  _ 

P 

In  the  incremental  form  this  condition  is: 

-AE - -A  p  —  (e  —  tn )  (5.22) 

P  P 

This  condition  is  imposed  via  the  penalty  method.  It  should  be  noted  that  there  are  available 
a  variety  of  possible  forms  of  the  penalty  terms,  depending  on  the  form  of  the  test  term 
applied  to  condition  (5.22).  One  possibility,  which  appears  to  be  the  most  natural  and  yields 
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a  symmetric  contribution  to  the  stiffness  matrix,  is  obtained  by  testing  equation  (5.22)  with 
its  own  variation: 


1 

e 


-A E  -  -A p 

P  P 


where  iq,)  denotes  a  selected  component  of  a  test  vector:  tqj)  for  density  and  v^\i)  (or  energy. 

This  approach  leads  to  two  penalty  conditions  affecting  both  the  continuity  and  energy 
equations.  Therefore  it  is  not  in  agreement  with  the  physical  situation  because,  while  the 
solid  wall  can  supply  heat  to  maintain  a  prescribed  termperature.  it  cannot  supply  mass  for 
this  purpose.  Due  to  this  reason,  another  form  of  the  penalty  term  should  be  used,  which 
enforces  a  prescribed  temperature  by  altering  the  energy  equation  only: 


1 

t 


(e-e") 


V(M) 


The  corresponding  terms  in  the  stiffness  matrix  and  the  right-hand  side  are: 

K,j  =  /  kVrVjds 

J  dQw 

Ri  =  I  r'ijds 

JdOw 


(5.23) 


where  the  kernel  matrices  k  and  r  are  defined  (in  two  dimensions)  as: 


k  = 


€ 


0  0  0  0 

0  0  0  0 

0  0  0  0 

E 

- 0  0  1 

p 


r  =  -(m,n ,) 
e 

Again,  nodewise  enforcement  of  these  conditions  can  be  obtained  by  replac  ing  sluqx1  (unc¬ 
tions  with  Dirac  delta  functions. 

For  the  regularized  problem  some  additional  artificial  terms  (fluxes)  occur  on  the  solid 
wall  due  to  second-order  terms  and  explicit  artificial  dissipation.  These  fluxes  are  forced  to 
be  zero  by  means  of  natural  boundary  conditions,  in  the  same  manner  as  for  the  supersonic 
outflow. 


0 

0 

0 

p(c-en) 


41 


5.5.6  No-Flow 


The  basic  condition  of  the  no-flow  boundary  is  that  the  normal  velocity  is  zero  or.  equiva¬ 
lently,  that  the  normal  momentum  is  zero: 

m,n,  -  0 


In  the  incremental  form  this  becomes: 


A  m,n,  =  — 


This  condition  is  easily  enforced  by  the  penalty  function,  with  the  additional  term  in  the 
variational  formulation: 


whe  re  {'(1+,)  is  the  test  function  corresponding  to  momentum  m,.  The  resulting  stillness  ma¬ 
trices  and  right-hand  sides  are  ot  the  standard  form  (5.23),  with  kernals  (in  two  dimensions): 

0  0  0  0' 

0  Mini  n\U2  0 

0  n2«l  0 

0  0  0  0 

1 

r  =  — 

£ 

For  the  viscous  flow  the  additional  conditions  on  the  no-flow  boundary  are  the  two  natural 
boundary  conditions: 


0 

Hr 

n-2 

0 


—  0 

<7n  =  0 

where  uns  is  the  skin  friction  on  the  boundary  and  qn  is  the  normal  heat  flux.  Tin?  appli¬ 
cation  of  these  natural  boundary  conditions  follows  the  procedure  discussed  in  preredinii 
subsections. 

Similarly,  as  on  the  solid  wall,  all  the  artificial  fluxes  are  forced  to  be  zero  on  the  no-tlow 
boundary. 
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5.6  Implicit/Explicit  Procedures 

The  basic  idea  of  implicit/explicit  algorithms  is  simple:  combine  the  two  methods  to  take  ad¬ 
vantage  of  superior  features  of  each  of  them.  The  major  advantage  of  the  explicit  method  is 
that  element  computations  are  relatively  cheap  and  simple.  Unfortunately  this  method  suf¬ 
fers  from  stability  limitations  of  the  time  step,  which  in  some  problems  leads  to  prohibitively 
large  numbers  of  time  steps. 

The  implicit  algorithm  is  unconditionally  stable  and  thus  allows  for  an  application  of 
larger  time  steps  than  explicit  method  (but  there  are  still  some  other  limitations).  More¬ 
over,  due  to  the  existence  of  implicit  boundary  terms  it  allows  for  an  easy  and  straightforward 
control  of  natural  boundary  conditions,  in  particular  the  viscous  fluxes.  An  additional  advan¬ 
tage  is  that  with  larger  time  steps  no  explicit  artificial  dissipation  is  required,  which  is  very 
important  in  the  calculation  of  boundary  fluxes,  in  particular  wall  heating  rates.  The  ma  jor 
disadvantage  is  much  higher  cost  of  element  operations  and  more  complex  and  expensive 
solution  of  the  resulting  system  of  equations. 

In  this  section  the  formulation  and  numerical  implementation  of  an  adaptive  implicit /ex¬ 
plicit  algorithm  will  be  presented.  The  algorithm  will  be  based  on  the  general  family  of 
Taylor-Galerkin  methods  discussed  in  the  previous  sections. 

5.6.1  Formulation  of  Implicit/Explicit  Schemes 

Assume  for  now  that  subregions  to  be  treated  with  implicit  and  explicit  methods  are  defined. 
Then  the  question  arises  as  to  what  is  a  proper  procedure  for  combining  implicit  /explicit 
computations.  "Proper"  here  means  that  it  maintains  conservative  properties  of  the  algo¬ 
rithm  and  the  required  order  of  approximation. 

Some  of  the  possible  procedures  are  examined  below: 

Procedure  I 

The  first  possible  approach,  applied  for  example  by  Hassan.  Morgan  and  Peraire  [17]  is  based 
on  the  following  two  steps: 

1.  Perform  the  explicit  step  computations  on  all  nodes  in  the  mesh. 

2.  Perform  the  implicit  computations  in  the  subregions,  where  the  stability  criterion  for 
the  explicit  scheme:  C  <  1,  is  violated.  The  solution  in  remaining  nodes  is  "frozen" 
at  this  step. 
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This  simple  procedure  has  one  basic  disadvantage:  it  appears  to  be  nonconservative  and 
it  may  disturb  the  regularity  of  the  solution  in  the  transition  zone.  This  is  caused  by  the 
fact  that  during  Step  2  the  "frozen"  explicit  nodes  impose  the  actual  Dirichlet  boundary 
condition  on  the  edge  of  the  implicit  zone.  Prescribing  Dirichlet  boundary  condition  for 
the  Xavier- Stokes  equation  means  that  there  must  exist  an  (external)  source  of  fluxes  to 
support  the  prescribed  state  of  the  solution.  Since  no  such  external  source  exists  within  the 
computational  domain,  the  solution  will  not  be  conservative  across  the  implicit /explicit  line. 
Due  to  enforcement  of  the  Dirichlet  boundary  conditions  the  solution  may  exhibit  a  "ramp ' 
or  "kink"  along  this  line. 

Procedure  II 

The  second  procedure  considered  is  based  on  a  division  of  the  computational  domain  H  into 
two  subregions  Q(/)  and  (see  Fig.  IT),  such  that: 

o(£)  n  f>(/)  =  0  .  P.{E]  U  P[I)  =  n 

It  is  convenient  to  assume  that  the  interface  between  the  two  regions  coincides  with  the 
element  boundaries. 

It  can  be  easily  observed  that  the  differential  equations  to  be  solved  on  the  two  subregions 
are  different  due  to  different  implicitness  parameters  applied  in  each  zone:  a*  .  dl/).  * !/*  in 
the  implicit  zone  and  o(El.  3{E) .  ~3E)  in  the  explicit  zone  (in  fact  the  latter  are  zero!.  There¬ 
fore.  the  variational  formulation  (5.14),  based  on  the  assumption  of  constant  implicitness 
parameters,  cannot  be  applied  to  the  domain  Q.  Instead,  it  can  be  applied  separately  to 
each  subdomain  with  additional  continuation  conditions  across  the  interface.  These  condi¬ 
tions  represent  continuity  of  the  solution  and  satisfaction  of  the  conservation  laws  across  the 
interface  and  are  of  the  form: 


=  u*7' 

=  F[n 

,{E) 

'  *  l 

=  A\n 

jp(E)V 

=  FU)V 

n 

where  index  n  refers  to  the  outward  normal  for  the  corresonding  region  i  n,/:l  =  The 

continuity  requirement  pertains  also  to  the  test  function,  so  that.  v(E]  =  u(/)  =  v.  If  the 
variational  statement  is  formulated  for  this  problem,  then  in  addition  to  interior  integrals  for 
each  subdomain  and  regular  boundary  integrals,  jump  integrals  across  the  interface  show  up 
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in  the  formulation.  These  additional  interface  terms  on  the  right-hand  side  are  of  the  form: 


-At  [/’J,/)n  4-  F{E)n]  ■  v  4-  At  \F[I)Vn  4-  F[E)Vn]  ■  v 


Atn 


~  [(l  -2a<'>)  A^F\']n  +  (1-2 a<E>)  A{E) F^)n\  •  vds 


It  can  be  easily  shown  that  the  first  two  terms  are  zero  because  of  the  interface  conditions 
(5.25).  On  the  left-hand  side  of  the  variational  formulation  additional  interface  terms  are  of 
the  form: 


-ra^  A^Au^n^ 


•  v 


-  At  jf  'I  (Rl'l  n  +  P[n  ■  A.Wn!0)  +  1|£|  (Rifa»«/>nP  +  PlfUU'E'..!E|)j  .  „ 


-  ~  [(,  _  3.(0)  +  (,_  2a!*)) 


(£)n(£) 


vds 


or,  after  reinterpretation  of  the  linearized  terms: 

At  +al£>^Fi,£l]  ■  v  -  AT  [-j{n  AF[I)V  -  AF[E]V]  v 

-  42  [u  _  iatn)3mjmAfin  +  (1  _  A<f' ■  Ws 

hi  the  above  formula,  natural  interface  conditions  should  be  imposed  on  the  viscous,  second- 
order  and  artificial  dissipation  terms.  Unfortunately,  basic  conditions  (5.25)  are  not  directly 
applicable  to  the  viscous  terms  because  of  different  coefficients  for  the  implicit  and  explicit 
terms.  Even  more  questionable  are  the  relevant  natural  interface  conditions  for  the  terms 
resulting  from  the  second-order  expansion  in  time.  At  present,  we  do  not  have  reasonable 
propositions  of  the  conditions  to  be  applied.  It  can  be  observed,  however,  that  all  the  terms 
considered  here  are  at  least  of  the  order  of  At*,  so  that  setting  them  all  to  zero  preserves 
the  first-order  time  accuracy  of  the  scheme.  This  is  the  approach  adopted  subsequently  in 
the  implementation  of  this  procedure. 


Procedure  III 

The  last  procedure  considered  here  is  based  on  a  generalization  of  the  weak  formulation, 
according  to  which  the  implicitness  parameters  are  not  constant,  but  are  continuous  functions 
of  the  position  x.  For  the  finite  element  computations  it  is  convenient  to  limit  the  choice  of 
these  parameters  to  the  finite  element  subspace,  so  that: 

.v 

a(j)  =  21 

i+ 1 
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and  the  same  holds  for  the  other  implicitness  parameters  j3, 7  and  8.  With  this  assumption 
the  weak  formulation  was  derived  for  the  domain  f l.  The  resulting  formula  is  of  the  form 
(5.14)  with  additional  terms  on  the  left-hand  side: 


L 


Ata  tA*  Au  ■  v 


+  A t  Aut]  +  7  ,P^Auj  ■  v 

-  All  [((l  _  oa)0).  ■  vdd 


and  on  the  right-hand  side: 


A  t1 


la  2 


■  vdfl 


From  these  equations,  the  formulas  for  additional  contributions  to  the  stiffness  matrix  and  the 
right-hand  side  can  be  easily  derived.  The  artificial  dissipation  contributions,  not  presented 
here,  are  handled  in  exactly  the  same  way  as  the  natural  viscosity.  Note  that  there  are  no 
additional  boundary  terms  resulting  from  the  variable  implicitness. 


The  above  approach  seems  to  be  the  most  general  and  clear,  with  no  ambiguities  con¬ 
cerning  the  interface  conditions.  To  make  the  practical  application  cheaper,  the  implicitness 
coefficients  are  held  constant  in  most  of  the  elements,  and  the  additional  terms  are  actually 
evaluated  only  in  the  transition  zones. 


5.6.2  Selection  of  Implicit  and  Explicit  Zones 

The  basic  criterion  for  selection  of  implicit  and  explicit  zones  is  simple:  for  a  given  time  step 
all  nodes  which  violate  the  stability  criterion  for  an  explicit  scheme  should  be  treated  with 
the  implicit  scheme.  According  to  this  criterion,  several  options  for  an  automatic  adaptive 
selection  of  implicit/explicit  zones  were  implemented: 

1.  User-prescribed  time  step  DT: 

Within  this  option  the  user  prescribes  the  time  step.  All  nodes  satisfying  stability 
criterion  for  the  explicit  scheme  (with  some  safety  factor)  are  explicit.  This  means 
that  all  the  elements  connected  to  these  nodes  are  treated  with  the  explicit  scheme. 
On  all  other  elements  the  implicit  scheme  is  applied. 

2.  Prescribed  maximum  CFL  number: 

In  this  option  the  user  prescribes  the  maximum  CFL  number  that  can  occur  for  ele¬ 
ments  in  Cl.  The  time  step  is  automatically  selected  as  the  maximum  step  satisfying 


46 


this  condition.  The  choice  of  maximum  CFL  number  may  be  suggested  by  the  time  ac¬ 
curacy  arguments  or  the  quality  of  results  (it  is  known  that  for  Taylor-Galerkin  scheme 
too  large  a  CFL  number  tends  to  smear  shocks — see  Section  6.3). 

3.  Prescribed  percentage  of  the  domain  is  implicit. 

In  this  version  the  user  specifies  the  fraction  of  the  domain  which  is  to  be  treated 
implicitly.  The  nodes  with  the  strongest  stability  limitation  (usually  the  smallest  ones) 
are  treated  implicitly,  the  others  are  explicit.  The  time  step  is  selected  to  guarantee 
stability  of  the  explicit  zone. 

4.  Minimization  of  the  cost  of  computations. 

In  this  option  the  time  step  and  the  implicit/explicit  subzones  are  selected  to  minimize 
the  cost  of  advancing  the  solution  in  time  (say  one  time  unit).  The  algorithm  is  based 
on  the  tact  that  for  an  increased  time  step  an  increasing  number  of  elements  must  be 
analyzed  with  the  (expansive)  implicit  algorithm.  The  typical  situation  is  presented  in 
Fig.  19.  which  shows  tor  different  time  steps  the  relative  number  of  nodes  that  must 
be  treated  with  the  implicit  scheme  (to  preserve  stability).  On  the  abscissa  the  Af/.-£ 
denotes  the  longest  time  step  allowable  for  the  fully  explicit  scheme  (with  some  safety 
factor).  Atfi  denotes  the  shortest  time  step  requiring  a  fully  implicit  procedure.  The 
number  of  implicit  nodes  increases  as  a  step  function  from  zero  for  At  <  At/.-/r  to  one 
lor  A/  >  tf  j.  Xow  assume  that  the  ratio  r  of  the  computational  cost  of  processing  one 
implicit  node  to  the  cost  of  processing  one  explicit  node  is  given  (including  the  time  of 
solving  the  system  of  equations).  This  ratio  is.  in  general,  a  function  of  the  number  of 
implicit  and  explicit  nodes.  Then  the  reduction  of  the  cost  of  advancing  the  solution 
in  time  with  the  implicit/explicit  scheme,  as  compared  to  the  fully  explicit  scheme,  is 
given  by  the  formula: 

R  =  („(£>  + 

Ai  '  ' 

Typical  plots  ot  the  {unction  R{t.)  are  presented  in  Fig.  19.  Shown  here  are  the  two 
cases: 


(a)  the  case  of  a  small  difference  between  fully  explicit  and  fully  implicit  time  steps — 
an  almost  uniform  mesh 

(1>)  the  case  of  a  large  difference  between  fully  explicit  and  fully  implicit  time  stops 

Before  discussing  these  plots  it  should  be  noted  that  certain  restrictions  on  the  length 
ot  the  time  step  should  be  applied,  for  example,  from  the  maximum  CFL  condition. 
Oth  erwise  the  cheapest  procedure  would  always  be  to  reach  the  final  time  with  one 
implicit  step. 
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From  the  plots  in  Fig.  19  the  following  observation  can  be  made:  for  an  almost  uniform 
mesh  the  mixed  implicit/explicit  procedure  does  not  provide  savings  of  the  computa¬ 
tional  cost — either  a  fully  implicit  or  fully  explicit  scheme  is  the  cheapest  depending  on 
the  time  step  restriction.  On  the  other  hand,  for  very  diverse  mesh  sizes  the  mixed  pro¬ 
cedure  provides  considerable  savings.  This  means  that  the  effectiveness  of  the  mixed 
implicit/explicit  scheme  will  be  the  best  for  large-scale  computations  with  adaptive 
meshes,  especially  problems  with  very  thin  boundary  layers.  Some  introductory  results 
are  presented  in  Section  6.3.  In  the  practical  implementation  of  this  method,  the  ap¬ 
proximation  to  the  function  R(At)  is  automatically  estimated  for  a  given  mesh.  Then, 
the  time  step  corresponding  to  the  smallest  /?(Af)  is  selected  automatically  (subject 
to  some  constraints,  in  particular  the  CFLmax  constraint). 

In  addition  to  the  above  few  criteria,  based  purely  on  a  stability  analysis,  some  other 
criteria  for  application  of  implicit  schemes  can  be  applied.  For  example,  within  boundary 
layers  the  implicit  scheme  may  be  preferred,  because  it  provides  faster  convergence  of  the 
boundary  fluxes  and  offers  direct  control  of  the  natural  boundary  conditions.  Many  of  these 
issues  are  yet  to  be  experimented  with. 

5.6.3  Computational  Procedure 

In  the  numerical  solution  of  the  Navier-Stokes  equations  the  implicit/explicit  procedure  is 
combined  with  an  adaptive  finite  element  mesh  refinement  scheme,  as  described  in  refer¬ 
ences  [28.29].  The  range  of  implicit/explicit  zones  is  redefined  after  every  mesh  refinement 
and — between  refinements — after  every  prescribed  number  of  steps.  At  every  time  step  the 
implicit/explicit  procedure  results  in  a  system  of  equations: 

(M  +  Kin)AU  =  R  to.26) 

where  the  stiffness  matrix  K  has  non-zero  entries  only  for  degrees  of  freedom  in  the  implicit 
zone  or  for  nodes  with  penalty-enforced  boundary  conditions.  The  solution  procedure  for  the 
above  system  differs  significantly  depending  whether  the  consistent  or  lumped  mass  matrix 
is  used. 

If  the  mass  matrix  is  lumped  (diagonal),  then  the  solution  is  performed  in  two  -tops: 

1.  calculate  the  solution  for  all  explicit  degrees  of  freedom  by  a  simple  one-pass  Jacobi 
procedure. 

2.  solve  the  linear  system  of  equations  for  the  implicit  degrees  of  freedom.  Currently  the 
two  solvers  are  available: 
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•  direct  frontal  solver  (for  testing  purposes), 

•  truncated-restarted  GMRES  method  with  block-Jacobi  preconditioning  and  an 
element-by-element  data  structure.  Other  preconditioners  are  currently  being 
tested. 

Note:  nodes  with  out-of-diagonal  penalty  terms,  like  for  example  the  ones  on  the  solid 
wall,  are  formally  classified  as  implicit  and  dealt  with  at  this  stage. 

If  the  mass  matrix  is  consistent,  then  in  spite  of  the  application  of  the  explicit  scheme,  the 
system  of  equations  is  implicit  (all  degrees  of  freedom  are  coupled).  Thus,  the  preconditioned 
GMRES  procedure  is  used  for  the  whole  system  (5.26). 

It  can  be  observed  that  the  implicit/explicit  procedure  is  computationally  much  more 
effective  if  the  lumped  mass  matrix  is  used.  The  explicit  degrees  of  freedom  are  eliminated 
at  very  low  cost  and  the  same  workspace  can  be  used  for  explicit  and  implicit  computations. 
The  consistent  mass  matrix  introduces  full  coupling  of  degrees  of  freedom  within  the  domain 
IT  Although  the  coupling  due  to  the  mass  matrix  is  rather  weak  and  would  require  just  a 
few  iterations,  full  consistency  of  the  solver  requires  iterating  in  the  explicit  zone  as  many 
times  as  in  the  implicit  zone.  This  considerably  increases  the  cost  of  solution  of  the  system 
of  equations,  especially  if  the  implicit  subzone  is  relatively  small. 

5.7  Computation  of  Aerothermal  Loads 

An  essential  task  in  the  numerical  analysis  of  hypersonic  fluid-structure  interactions  is  a 
precise  calculation  of  boundary  quantities,  in  particular  fluid  pressure,  skin  friction  and 
heating  rates.  These  quantities  contribute  to  the  mechanical  and  thermal  loads  exerted  by 
the  fluid  on  the  structure. 


5.7.1  Boundary  Fluxes 


Evaluation  of  aerothermal  loads  requires  the  computation  of  flux  across  the  boundary  f  ol 
the  fluid  domain  0.  namely: 

F„  =  F,n,  i  =  1.2  1 5.27 1 

where  F„  is  the  flux  vector  across  the  boundary,  n  is  the  outward  normal  to  the  boundary 
T  and  F,  are  flux  vectors  in  directions  of  the  global  coordinate  axes  ,r,  (see  Fig.  5).  Each  of 
these  vectors  can  be  presented  as  a  composition  of  inviscid  and  viscous  fluxes. 
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with  components  of  the  vectors  F,  and  F\-  given  in  equation  (5.4). 
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The  normal  flux  can  also  be  decomposed  into  inviscid  and  viscous  parts: 

XT'  —  TTC  _  J?V 
r  n  —  -C  n  n 

which,  after  application  of  equation  (5.27),  are  represented  as: 


(5.2S) 
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[  vktl£  -  qn 

where  vn  denotes  normal  outward  velocity  and  t)  are  components  of  the  boundary  traction 
due  to  viscous  forces. 

The  consecutive  components  of  the  normal  flux  vector  Fn  can.  in  general,  be  interpreted 
as  mass  flux,  .rj- momentum  flux.  x>-momentum  flux  and  total  energy  flux  across  the  bound¬ 
ary. 

On  a  solid  wall,  which  is  of  primary  interest  in  our  analysis,  the  components  of  the 
velocitv  vector  vanish,  and  the  total  boundary  flux  is: 


Fa 


(5.30) 
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In  this  case  the  total  normal  pressure  on  the  wall  p  and  the  value  of  the  skin  friction  r  are 
calculated  as  projections  of  the  xq  and  xq-momentum  fluxes  on  vectors  n  and  s.  which  are 
normal  and  tangent  to  the  boundary,  respectively: 

P  =  Fn{2)n\  4  Fn^n2  =  p  —  cr\n 


7  -  Fn[2)S\  4-  F ^(3)^2  -  -ffL 

Note  that  the  normal  pressure  on  the  wall  is.  in  general,  different  than  the  thermodynamic 
pressure.  It  can  be  verified  that  in  laminar,  fully  developed  boundary  layers  viscous  contri¬ 
butions  to  the  normal  pressure  cr^n  vanish,  but  they  can  be  significant  in  other  cases,  like 
the  shock-boundary  layer  interaction. 

Regarding  practical  evaluation  of  boundary  fluxes  in  the  finite  element  method,  there  are 
two  basic  approaches,  namely: 
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-  calculation  from  the  definition. 

-  consistent  finite  element  calculation. 

The  calculation  from  the  definition  is  straightforward  and  reduces  to  a  pointwise  application 
of  formula  (5.4).  The  necessary  gradients  of  components  of  velocity  and  the  temperature  are 
calculated  using  the  finite  element  shape  functions.  Since  these  gradients  are  discontinuous 
across  element  boundaries,  the  nodal  values  of  boundary  fluxes  can  be  obtained  by  averaging 
of  the  values  from  the  adjacent  element  sides  (or  by  other  postprocessing  procedures). 

The  calculation  of  boundary  fluxes  from  the  definition  has.  for  most  commonly  used  linear 
or  bilinear  elements,  poor  accuracy  (first  order  in  space).  More  precise  results,  corresponding 
approximately  to  the  application  of  second  order  schemes  in  the  finite  difference  method, 
can  be  obtained  from  a  consistent  finite  element  recovery.  This  method  will  be  presented  in 
the  next  section. 


5.7.2  Consistent  Computation  of  Boundary  Fluxes 


The  basis  for  a  consistent  calculation  of  boundary  fluxes  is  The  weak  formulation  of  the 
Xavier-Stokes  equations,  written  in  the  form: 

[  u  ■  vd9.  —  /  Ft  ■  v.,d9.  +  [  Fn  vds  =  0  V  v  £  V"  (5.311 

/  n  J  n  J'.'O 

where  V’  is  the  space  of  test  functions  (usually  a  subspace  of  H](Cl)). 

Suppose  that  the  domain  Q  was  discretized  with  finite  elements,  the  approximate  solution 
of  equation  (5.31 )  was  obtained  and  we  want  to  calculate  boundary  fluxes  Fn  across  a  certain 
section  of  the  boundary  F  f-  (see  Fig.  6). 

The  general  idea  is  to  reformulate  equation  (5.31)  in  the  form: 

f  Fn  ■  vdS  =  —  /  it  ■  vdP.  f  F iV  ,f/n  —  /  Fn  ■  vdS  (5.32) 

J  rF  Jci  Jo.  Ji‘0.-rr 

The  test  functions  in  this  equation  are  arbitrary  linear  combinations  of  nodal  shape'  functions 
'P/.  For  the  purpose  of  calculating  fluxes  on  F/r  it  suffices  to  restrict  the  domain  (.)  to  the 
subdomain  91.  which  is  the  set  of  elements  covered  by  shape  functions  associated  with  nodes 
on  the  boundary  I'/-  (hatched  areas  in  Fig.  21).  I  hen  the  test  functions  ronsideied  are  •  .1 
the  form: 

X  f 

v(x)  =  ^  V i(x)  (5.33) 

t= i 

Thus  equation  (5.32).  after  using  the  fact  that  v  =  0  for  the  boundary  F.  takes  the  form: 

/  Fn  ■  vdS'  =  —  fu  ■  vd9l  —  [_  Ft  ■  v,dQ  —  /  Fn  ■  vdS  (5.34) 

JvF  J n  J n  Jvo 
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Introducing  the  finite  element  approximation  of  fluxes  along  the  boundary  T p 


‘VF 

Fn(x )  =  FniWj(x) 

/  =  1 

and  requiring  that  equation  (5.34)  be  satisfied  for  every  test  function,  leads  to  the  following 
linear  system  of  equations: 

MljFnJ  =  Rj 

In  the  above.  M  is  the  "mass"  matrix  for  the  boundary  Ff. 

Mu=f  I^V^VjdS 
J  r> 

where  IAxA  represents  unit  matrix  of  size  4x4.  For  simplicity,  a  lumped  mass  matrix  can 
also  be  used.  The  right-hand  size  vector  is  calculated  according  to  the  formula: 

R.j  =  j -r  F,y]> j,tdP. 

-  /  Fn'ltjdS 
J  r0 

Note  that  it  is  essential  to  include  in  the  calculation  of  the  right-hand  sides  the  contribution 
of  boundary  terms  on  Fq.  If  these  terms  are  neglected,  then  incorrect  values  of  lluxes  are 
calculated  near  both  ends  of  Ff.  The  necessary  fluxes  F n  on  F0  are  calculated  in  the  same 
way  as  the  fluxes  F,  inside  0.  namely  using  definition  (5.27)  and  the  gradients  calculated 
from  finite  element  shape  functions. 

It  is  important  to  observe  that  the  application  of  artificial  dissipation  models  (like  Lapidus 
dissipation!  introduces  additional  terms  into  the  equation  (5.31).  These  terms  give  rise  to 
artificial  boundary  fluxes,  which  have  no  physical  interpretation.  Thus,  including  these  terms 
in  computations  yields  results  of  no  physical  meaning.  On  the  other  hand,  neglecting  these 
terms  in  heat  flux  computations  violates  satisfaction  of  the  weak  form  of  the  conservation 
laws  (5.31  I.  so  the  calculated  heat  fluxes  are  not  consistent  with  the  finite  element  solution. 
Therefore,  the  conclusion  is  that,  for  the  purpose  of  a  propci  evaluation  of  boundary  fluxes, 
the  artificial  dissipation  should  not  be  used  in  the  vicinity  of  a  solid  wall  ''at  least  in  the 
present  form  (5.15/  or  |5.!t«i). 

6  Selected  Validation  Examples 

In  order  to  better  understand  phenomena  of  hypersonic  interactions  and  t  o  validate  t  be  met  h- 
ods  and  software  developed  in  the  project,  several  test  problems  were  solved.  The  problems 
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selected  were  usually  rather  simple,  so  that  experimental  results,  analytical  solutions,  or 
solutions  obtained  by  other  researchers  could  be  identified  for  comparison.  Most  of  these 
validation  examples  were  presented  in  consecutive  progress  reports  and  will  not  be  repeated 
in  this  report.  However,  in  order  to  illustrate  the  correctness  of  the  proposed  formulation, 
methods  and  software,  some  of  these  results  will  be  briefly  summarized.  Moreover,  recent 
examples  of  implicit/explicit  flow  computations  will  be  presented  in  this  section. 

6.1  Validation  of  Thermo— Viscoplastic  Analysis 

The  first  simple  verification  test  of  the  viscoplastic  constitutive  equations  is  the  simple  tension 
test  at  elevated  temperatures,  performed  for  a  superalloy  B1900  -fHf  under  the  NASA  HOST 
program  [5.6].  As  indicated  in  Fig.  22.  at  1S00°F  the  elastic  section  of  a  stress-strain  curve 
is  very  short  and  the  difference  between  elastic  and  viscoplastic  solutions  is  very  significant. 
The  predictions  of  the  Bodner-Partom  unified  viscoplastic  theory  correspond  very  well  to 
experimental  results. 

To  verify  the  characteristics  of  the  thermo-viscoplastic  constitutive  models  under  more 
complex  load  histories,  the  stress-strain  response  of  a  bar  due  to  specified  cyclical  variations 
of  the  strain  rate  was  analyzed.  Figure  23  presents  a  statement  of  the  problem,  and  Fig.  24 
shows  the  stress-strain  response  of  the  bar  for  different  temperatures.  Significant  points  to 
note  are  that  with  increasing  temperature,  the  reduction  in  elastic  modulus  and  yield  stress 
are  correctly  represented  by  the  unified  viscoplastic  theory.  The  above  results  compare 
favorably  with  experimental  observations  for  nickel  superalloy  B1900  -t-Hf. 

6.2  Validation  of  Large  Displacement  Formulation 

The  large  displacement  extension  was  introduced  into  the  theory  and  software  in  order  to 
properly  capture  the  beam-column  effect  due  to  localized  heating  of  hypersonic  structures. 
To  validate  this  modification  and  to  verify  the  adequacy  of  the  finite  element  model,  an 
elastic  beam-column  problem  was  analyzed.  The  problem  statement  is  presented  in  Fig.  2b. 
A  simply  supported  beam  is  subjected  to  a  sudden  increase  in  surface  heating  on  its  upper 
surface.  The  lower  surface  of  the  beam  is  adiabatic.  The  heating  which  is  uniform  over 
the  upper  surface'  of  the  beam  induces  a  time-varying  temperature  variation  through  the 
beam  depth.  For  supports  fixed  against  horizontal  displacement,  the  beam  experiences  an 
increasing  axial  compressive  force.  At  the  same  time,  the  temperature  gradient  through  the 
depth  induces  a  thermal  moment,  and  so  the  beam  displaces  transversly  as  a  beam-column. 
1  he  elastic  beam-column  problem  has  a  closed-form  analytical  solution  that  can  be  used  to 
validate  the  finite  element  analysis. 


The  particular  data  used  in  this  analysis  were 


Q  =  6.0  BTU/s-in' 
p  =  0.2S3  1  b  m  /  i  n  3 
c  =  0.20  BTU/lbm  -  R 
k  =  4.16  •  10-3  BTL’/in-sec  —  R 
E  =  1.47  •  105MPa 
v  =  0.1 

a  =  0.S61 1  •  10~5  1  /R 

where  0  is  the  prescribed  heat  flux,  p  is  the  density  of  a  material,  c  is  the  thermal  capacitance. 
k  is  the  thermal  conductivity.  E  is  the  elastic  modulus,  v  is  the  Poisson's  ratio,  and  a  is  the 
thermal  expansion  coefficient. 

The  above  problem  was  solved  using  three  different  finite  element  structural  models. 
These  models  and  deformed  configurations  of  the  beam  at  t  =  1.0s  are  shown  in  Fig.  26. 

A  comparison  of  the  finite  element  solutions  with  the  beam-column  theory  is  presented 
in  Fig.  27.  Axial  force  versus  transverse  deflection  predictions  for  the  three  finite  element 
solutions  are  compared  with  the  theory.  The  results  for  the  most  refined  mesh  validate 
the  finite  element  large  displacement  theory.  The  results  also  show  that  the  cruder  meshes 
predict  too  large  an  axial  force  with  forces  becoming  higher  than  the  Euler  buckling  load  as 
the  deflection  increases  (errors  in  deflection  of  the  beam  were  relatively  small-order  of  a  few 
percent). 

6.3  Validation  of  Damage  Modeling 

To  gain  insight  into  behavior  of  a  viscoplastic  structure  with  damage  represented  in  the 
constitutive  equations  a  simple  model  was  analyzed  (Fig.  2$).  The  model  was  subjected  to 
ten  consecutive  tension-compression  cycles  of  strain  for  three  specified  temperature  levels. 
The  properties  of  the  material  correspond  to  the  superalloy  B 19000  4-Hf  as  presented  in 
Refs.  [5]  and  [6j. 

Figure  29  shows  st  r  —  t  ra  in  curves  with  damage  evolution  for  three  temperature  levels. 
At  the  lower  temperature  levels  and  tor  the  prescribed  range  ot  strain  rates,  i  here  is  little 
apparent  effect  of  the  damage,  but  at  the  higher  temperature  (1900°  C)  there  is  gradual 
reduction  in  the  yield  stress  with  increasing  strain  cycles. 
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6.4  Implicit/Explicit  Analysis  of  Compressible  Flow 


Several  examples  of  compressible  flow  were  solved  using  the  implicit/explicit  algorithm  de¬ 
veloped  in  Section  5.  The  purpose  of  these  analyses  was  to  verify  the  performance  of  the 
scheme  as  well  as  to  compose  explicit  and  implicit  versions  of  the  algorithm  and  to  identify 
possible  indicators  for  the  selection  of  implicit  and  explicit  subzones. 

Several  different  examples  were  solved  with  various  selections  of  parameters.  Some  of 
these  examples  will  be  presented  in  this  report. 

6.4.1  A  Shock  Tube  Problem 

A  shock  tube  problem  is  a  good  verification  example  for  transient  algorithms.  A  tube  of 
length  50  is  divided  by  a  diaphragm  into  two  sections,  with  fluid  at  rest  at  the  following 
values  of  density  and  pressure: 

pi  —  10.0  pi  =  7.1429 

P2  =  1-0  P2  =  0.71429 

After  the  diaphragm  is  removed,  the  transient  flow  developes  with  such  phenomena  as  a 
shock,  a  contact  region  and  an  expansion  region. 

The  above  problem  was  solved  with  different  combinations  of  parameters  and  the  results 
were  compared  with  the  analytical  solution.  The  comparison  of  the  density  distribution 
after  100  time  steps  with  a  CFL  =  0.25  (total  time  13.16S)  for  selected  combinations  of 
input  parameters  are  presented  in  Fig.  30.  The  general  conclusion  of  these  (and  other) 
results  are- 

•  with  small  time  steps  (CFL  <  1)  the  explicit  and  implicit  algorithms  give  the  same 
results. 

•  for  small  time  steps  (CFL  <  1)  artificial  dissipation  is  necessary  to  suppress  oscillations 
of  the  solution. 

•  larger  time  steps  (CFL  >  5)  combined  with  artificial  dissipation  tend  to  excessively 
smear  discontinuities,  and 

•  tor  larger  time  steps  (CFL  >  5)  no  explicit  artificial  dissipation  is  necessary  to  suppress 
excessive  oscillations  of  the  solution  (except,  for  the  diaphragm  zone  in  this  example). 
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6.4.2  Inviscid  Compression  Corner 


Another  typical  test  problem  is  a  steady  state  solution  of  a  compression  corner.  The  par¬ 
ticular  flow  considered  is  that  of  Mach  3  and  the  ramp  angle  16°.  The  problem  was  solved 
using  adaptive  refinement  techniques  with  fully  implicit  and  explicit  versions  of  the  algo¬ 
rithm.  The  results  obtained  with  both  these  versions  at  the  same  time  step  (CFL  =  0.5)  are 
practically  the  same.  The  refined  mesh  and  density  contours  are  presented  in  Fig.  31.  It  ran 
be  observed  that  the  algorithm  developed  in  this  work  floes  not  generate  any  artificial  shock 
reflections  on  the  outflow  boundary  (typical  for  many  implementations  of  Taylor-l  lalerkin 
or  Lax-WendrofT  algorithms). 

The  next  figure  presented  is  the  solution  obtained  with  implicit  algorithm  for  CFL  = 
10  with  no  explicit  artificial  dissipation.  It  can  be  observed  that  although  some  oscillation 
of  the  solution  can  be  observed,  the  shock  is  captured  within  two  elements.  This  result 
confirms  the  fact  that,  for  larger  time  steps,  the  algorithmic  diffusion  due  to  second-order 
terms  is  sufficient  to  suppress  spurious  oscillations  and  there  is  no  need  for  explicit  artificial 
dissipation. 

6.4.3  Boundary  Layer  Problem 

To  analyze  the  behavior  of  implicit  and  explicit  schemes  in  the  boundary  layer  the  supersonic 
boundary  layer  problem  shown  in  Fig.  33  was  solved.  The  problem  consists  of  specjfvim> 
inflow  profiles  of  conservation  variables  and  computing  profiles  at  downstream  locations  by 
solving  the  Xavier-Stokes  equations.  The  data  chosen  is  representative  of  a  supersonic  How. 
and  the  temperature  profiles  iFig.  34)  are  reprsentative  of  flows  that  induce  aerodynamic 
heating.  As  indicated  by  the  computations  presented  in  previous  reports,  the  artificial  dissi¬ 
pation  has  to  be  "turned  off"  in  the  boundary  layer  in  order  for  hearing  rates  to  be  predicted 
properl}-.  Recently  this  problem  was  solved  with  the  implicit  algorithm  at  different  time  steps. 
The  general  conclusion  is  that  within  the  boundary  layer  the  heating  rates  with  largo  tilin' 
steps  (CFL  numbers  up  to  500)  are  virtually  the  same  as  for  small  time  steps. 

6.4.4  Flat  Plate  Problem 

1  lie  problems  presented  so  far  were  solved  using  the  fully  explicit  ,,(■  fully  jmpiiejt  v.-jA<>n  of 
the  algorithm.  As  a  test  problem  for  the  implicit/explicit  procedures,  a  classical  flat  plate 
problem  solved  originally  by  Carter  was  analyzed.  The  problem  statement  for  this  problem 
is  presented  in  Fig.  35.  the  density  contours  corresponding  to  the  stead}-  state  solution 
arc  presented  in  Fig.  36.  and  the  nondimensional  heating  rates  on  the  wall  calculated  as 
'/"  =  '//fix.t'L  presented  in  Fig.  37. 


In  the  solution  of  the  problem  the  implicit/explicit  algorithm  was  used.  The  selection  of 
implicit  and  explicit  zones  was  based  on  the  minimization  of  the  cost  of  computations  with 
the  limitation  on  the  time  step:  CFL  <  5.0.  The  safety  coefficient  for  the  explicit  method  was 
equal  0.7.  For  the  mesh  presented  in  Fig.  36.  the  automatically  selected  nondimensional  time 
step  was  equal  4.17073.  which  corresponds  to  CFL  =  1.95.  The  number  of  explicit  elements 
was  equal  1327  and  the  number  of  implicit  elements  was  equal  90.  For  these  parameters 
the  cost  of  advancing  the  solution  in  time  was  1.6  times  cheaper  than  with  the  fully  explicit 
scheme  and  about  ten  times  cheaper  than  the  fully  implicit  scheme. 

The  above  analysis,  although  performed  with  a  very  conservative  restriction  on  the  time 
step,  confirms  the  effectiveness  of  adaptive  implicit/e.xplicit  scheme  in  the  solution  of  N'avier- 
Stokes  problems.  Further  experiments  are  necessary  to  refine  the  criteria  for  implicit /explicit 
calculations  and  further  increase  the  effectiveness  of  this  approach. 

The  fiat  plate  example  was  solved  primarily  to  test  the  implicit/e.xplicit  procedure  with 
no  special  emphasis  placed  on  obtaining  the  best  quality  results.  It  can  be  seen,  however, 
that  for  the  present  mesh  the  general  flow  phenomena  and  boundary  conditions  are  properly 
resolved.  In  particular,  the  fine  mesh  near  the  tip  of  the  plate  satisfactorily  captures  the 
singular  character  of  the  solution.  However,  a  few  critical  remarks  about  these  results  are 
appropriate  here: 

•  The  automatic  mesh  refinement  procedure  based  on  the  density  error  estimate  does  not 
sufficiently  refine  the  boundary  layer  region.  More  general  error  estimators,  including 
velocity  gradients,  have  recently  been  implemented  to  resolve  this  problem. 

»  The  heating  rate  plot  shows  some  disturbance  in  the  smoothness  of  the  graph  at 
points  of  changing  element  size.  This  is  due  to  the  application  of  artificial  dissipa¬ 
tion  throughout  the  whole  computational  domain.  Since  the  significance  of  the  artifi¬ 
cial  dissipation — as  compared  with  natural  dissipation— depends  on  the  mesh  size,  the 
rapidly  varying  element  size  results  in  a  disturbance  of  the  heating  rates.  As  indicated 
in  the  previous  section,  a  proper  computation  of  the  heating  requires  the  "turning  of" 
the  artificial  dissipation  in  the  boundary  layer  (or  extremely  small  elements). 

It  should  be  emphasized  at  this  point,  from  an  examination  of  the  right-hand  side  of 
the  variational  formulation  (5.14).  that  a  steady  state  solution  for  a  given  turn  >l<p  dots 
not  depend  on  the  implicitness  or  explicitness  of  the  solution  algorithm  (for  unconditionally 
stable  versions  of  the  implicit  scheme).  The  only  exceptions  are  boundaries  with  prescribed 
natural  boundary  conditions,  where  the  implicit  algorithm  provides  direct  control  of  the 
fluxes. 
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7  Thermo-Viscoplastic  Analysis  of  a  Convectively 
Cooled  Structure 

A  more  realistic  two-dimensional  model  of  a  convectively  cooled  hypersonic  structure 
is  shown  in  Fig.  3S.  The  model  represents  a  segment  of  a  convectively  cooled  structure 
such  as  a  wall  of  a  scramjet  engine  fuel  injection  strut.  The  finite  element  thermal  model 
(Fig.  3$)  includes:  (1)  conduction  heat  transfer  in  the  aerodynamic  skin,  heat  exchanger 
fins  and  a  primary  structure.  (2)  convective  heat  transfer  between  the  walls  of  the  coolant 
passage  and  coolant.  (3)  mass  transport  convection  in  the  coolant  which  has  an  unknown 
bulk  temperature,  and  (4)  surface  radiation  on  the  aerodynamic  skin.  The  aerodynamic  skin 
is  uniformly  convectively  heated  over  its  length,  and  superimposed  on  the  uniform  heating 
is  a  local,  intense  heating  simulating  a  transient  shock.  The  transient  heating  is  induced  by 
the  time-dependent  convection  coefficient  shown  in  Fig.  3S. 

In  the  plane  strain  structural  finite  element  model,  the  primary  structure  and  aerody¬ 
namic  skin  have  unit  thickness,  but  the  heat  exchanger  fins  in  the  coolant  passage  are  ap¬ 
proximately  represented  by  a  single  fin  with  the  total  thickness  of  0.060.  The  wall  segment 
has  fixed  displacement  boundary  conditions  at  the  left  and  right  ends.  The  aerodynamic 
skin,  exchanger  fins  and  primary  structure  are  also  loaded  by  internal  pressure  (7.0  MPa  = 
1000  psi)  in  the  coolant  passage  which  is  relatively  large  compared  to  external  aerodynamic 
pressure  (0.3-5  MPa  =  50  psi). 

The  model  was  analyzed  first  for  steady  aerodynamic  heating  and  internal  pressure  op¬ 
erating  conditions.  Then  the  sudden  localized  aerodynamic  heating  was  applied.  Thus, 
the  steady  temperatures  and  stresses  serve  as  initial  conditions  for  the  transient  thermo¬ 
viscoplastic  ana.  sis. 

The  viscoplastic  solution  was  computed  with  the  variable  time  step  algorithm.  A  total 
of  265  steps  were  required  to  compute  the  response  for  a  total  time  duration  of  1.2  s. 

The  thermal  response  of  the  wall  segment  is  shown  in  Fig.  39.  Figure  39b  shows  the 
time  history  of  the  temperature  at  a  point  on  the  aerodynamic  skin  directly  under  the 
transient  heating:  Fig.  39a  shows  contours  of  temperatures  at  t  —  0.5s  when  temperatures 
are  maximum.  The  temperature  history  is  qualitatively  similar  to  the  results  obtained  in  the 
one-dimensional  model  showing  the  rapid  skin  temperature  rise  and  fall  with  the  simulated 
shock  heating.  The  temperature  contours  show  relatively  steep  thermal  gradients  induced 
by  the  local  heating,  and  that  the  coolant  substantially  limits  the  extent  of  the  induced  high 
temperatures.  Thus  the  high  temperatures  are  confined  to  the  aerodynamic  skin,  and  the 
primary  structure  experiences  only  small  temperature  changes.  The  temperature  gradients 
in  the  skin  particularly  at  the  coolant-skin  interface  are  not  predicted  with  high  accuracy 
because  of  the  engineering  model  of  the  coolant  heat  transfer.  However,  local  temperature 


58 


levels  in  the  skin  are  reasonably  accurate  since  the  net  energy  transfer  to  the  coolant  is 
modeled  satisfactorily. 

Histories  of  the  horizontal  stress  component  ox  at  two  points  through  the  skin  thickness 
are  shown  in  Fig.  40.  The  stress  histories  follow  the  temperature,  and  under  the  intense  local 
heating  stresses  are  very  similar  to  the  results  obtained  from  the  one-dimensional  model. 
At  this  location  the  skin  yields  through  most  of  its  thickness.  After  the  heating  ceases  there 
is  a  rapid  decay  of  stress,  and  under  the  intense  local  heating  there  are  residual  tensile 
stresses.  Figure  41  shows  the  time  history  of  the  vertical  component  of  stress  a.,  in  the 
heat  exchanger  fins.  High  tensile  stresses  are  induced  with  significant  local  yielding  which 
is  followed  by  residual  compressive  stresses.  The  crude  finite  element  model  of  the  fins  only 
approximates  these  stresses,  but  the  high  tensile  stresses  can  potentially  cause  a  bond  failure 
at  the  heat  exchanger/skin  joint.  The  consequences  of  such  a  failure  are  studied  in  the  next 
example. 

The  viscoplastic  stress  histories  are  compared  with  stresses  predicted  assuming  an  elas¬ 
tic  behavior  in  Figs.  42  and  43.  The  elastic  computations  were  made  with  temperature- 
dependent  elastic  properties.  The  loss  of  stiffness  due  to  the  elevated  temperatures  for 
04/  <  t.  <  0.5s  accounts  for  the  stress  "rolling-over'’  with  an  appearance  of  yielding.  Gen¬ 
erally,  the  elastic  stresses  are  too  high  and.  of  course,  since  yielding  does  not  occur  residual 
stresses  are  not  predicted.  Maximum  deformations  (not  shown)  predicted  by  the  elastic 
analysis  and  viscoplastic  analysis  are  about  the  same,  approximately  0.001  in. 

The  two-dimensional  character  of  the  stress  components  ax  and  <7y  are  shown  in  Figs  44 
and  45.  respectively.  The  stress  distributions  are  shown  at  three  times  t  =  0.05  s.  0.5s.  1.2s 
in  the  response.  The  stresses  at  t  =  1.2s  are  the  residual  stresses.  Figure  46  shows  contours 
of  the  principal  plastic  strains  at  t  =  0.5s.  The  contours  show  the  relatively  localized  nature 
of  the  high  stress  gradients  as  well  as  the  tensile  and  compressive  regions.  The  significant 
residual  stresses  suggest  the  possibility  of  cumulative  damage  under  repeated  load  cycles. 
This  possibility  needs  further  investigation. 

Convectively  Cooled  Structure  With  Damage 

The  preceding  analysis  showed  high  tensile  stresses  in  the  heat  exchanger  fins  and  suggested 
the  possibility  of  an  aerodynamical!}’  skin/heat  exchanger  bond  failure.  To  simulate  such 
a  possibility,  the  analysis  was  repeated  but  with  the  fins  given  greatly  reduced  material 
properties  for  a  width  of  0.30  in.  under  the  high  localized  heating.  !  wo  results  of  this 
analvsis  are  shown  in  Figs  47  and  48. 

Since  the  fins  can  no  longer  support  tensile  stress,  the  internal  coolant  pressure  and  the 
beam-column  effect  causes  substantial  local  deformation  of  the  aerodynamic  skin,  figure  4  i 
shows  the  deformed  structure  at  t  =  0.5s.  Compared  are  the  solutions  obtained  with  small 
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deformation  and  large  deformation  theory.  Deformations  are  shown  to  scale,  and  the  figure 
shows  that  the  large  displacement  formulation  predicts  larger  displacements.  This  compar¬ 
ison  is  made  quantatively  in  Fig.  4S.  Viscoplastic  displacement  histories  are  compared  on 
the  aerodynamic  skin  for:  (1)  fins  without  debonding,  (2)  fins  with  debonding  as  predicted 
by  the  infinitesimal  formulation,  and  (3)  fins  with  debonding  as  predicted  by  the  large  dis¬ 
placement  formulation.  The  results  clearly  show  that  the  large  displacement  formulation  is 
required. 


8  Conclusions 

In  the  course  of  this  project  considerable  progress  was  made  towards  the  understanding 
of  complex  phenomena  of  flow-thermal  and  structural  interactions  in  hypersonic  flight  and 
towards  reliable  methods  of  numerical  modeling  of  these  phenomena.  The  result  of  the 
project  is  a  unique  computational  capability  for  modeling  hypersonic  flow,  prediction  of 
aerothermal  loads  and  a  thermo-viscoplastic  analysis  of  the  structural  response. 

The  numerical  analysis  of  all  three  components  of  the  interaction  problem  (flow,  thermal, 
and  structural)  is  based  on  the  finite  element  method,  thus  providing  full  compatibility 
and  coupling  of  these  analyses.  To  guarantee  efficient  computations  and  reliable  results, 
the  software  is  equipped  with  adaptive  mesh  refinement,  adaptive  time  step  selection,  and 
adaptive  implicit/explicit  procedures,  based  on  relevant  error  and  cost  estimates. 

An  analysis  of  a  two-dimensional  model  of  a  realistic  convectively  cooled  structure  pro¬ 
vided  more  detailed  understanding  of  the  thermal-structural  behavior.  The  coolant  flow 
dominates  the  thermal  response  providing  a  relatively  short  thermal  transient.  I'nder  in¬ 
tense  heating,  significant  local  plasticity  occurs  in  the  aerodynamic  skin,  but  the  primary 
structure  remains  undamaged.  Heat  exchanger  fins  experience  high  tensile  stresses  and  both 
the  aerodynamic  skin  and  heat  exchanger  fins  have  significant  residual  stresses.  An  analysis 
of  an  aerodynamic  skin/heat  exchanger  fin  bond  failure  showed  that  the  aerodynamic  skin 
would  experience  significant  local  plastic  deformation  due  to  combined  coolant  pressure  and 
beam-column  effect  in  the  skin.  The  local  deformation  is  pronounced  enough  to  disturb  the 
external  aerodynamic  flow  and  introduce  a  complex  fiow-therma'-structural  interaction. 

The  analysis  of  convectively  cooled  hypersonic  structurs  under  intense  local  heating  shows 
the  important  role  that  such  analyses  can  make  in  understanding  complex  transient  inelastic 
structural  behavior  at  elevated  temperatures. 

Further  progress  in  the  realistic  modeling  of  hypersonic  flow-thermal  and  structural  in¬ 
teractions  will  involve  both  theoretical,  numerical,  and  experimental  analyses.  Some  of  the 
issues  to  be  resolved  are: 
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1.  An  experimental  verification  of  thermo-viscoplastic  structural  analyses.  Although  the 
viscoplastic  Bodner-Partom  equations  compare  very  favorably  with  experiments  for 
simple  stress  states,  extensive  experimental  verification  is  necessary  to  confirm  their 
adequacy  in  complex  stress  states  due  to  localized  heating. 

2.  Further  refinement  of  computational  techniques  towards  maximum  efficiency  and  re¬ 
liability  is  needed.  The  hypersonic  interaction  problems  are — due  to  combined  tluid- 
thermal-structural  analyses  and  due  to  their  inherent  transient  nature — computation¬ 
ally  very  expensive.  Further  progress  toward  computational  efficiency  is  necessayr  to 
make  realistic  three-dimensional  analyses  feasible. 

3.  Fast  and  reliable  techniques  for  coupling  of  structural  deformations  with  the  flow  anal¬ 
ysis  are  needed.  The  issues  include  moving  boundaries,  transfer  of  aerothermal  loads, 
matching  of  time  scales,  etc. 

4.  Constitutive  models  for  describing  high-temperature  behavior  of  nonmetallic  materials 
should  be  developed  and  implemented.  For  metallic  materials,  constitutive  relations 
valid  in  temperatures  near  the  melting  point  should  be  developed. 

■5.  Models  and  computational  techniques  for  surface  reactions,  burn-through  and  other 
failure  phenomena  need  to  be  developed  and  implemented. 

6.  Extensive  numerical  analyses  and  experimental  verifications  are  necessary  to  enhance 
the  understanding  of  phenomena  of  hypersonic  interactions  and  verify  the  reliability 
of  theoretical  models  and  computational  techniques. 

Resolution  of  several  of  these  issues,  in  particular  an  experimental  verification  of  thermo- 
viscoplastic  numerical  analyses,  were  planned  for  the  third  year  of  this  effort. 
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Figure  4:  Engineering  model  of  coolant  passage  heat  transfer. 


Figure  5:  Finite  element  model  of  coolant  passage. 
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Figure  7:  Large  deformation  of  a  deformable  continuum. 
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Figure  10:  Stability  analysis  for  fully  explicit  second-order  scheme. 


Figure  11:  Stability  analysis  tor  weakly  implicit  scheme:  a  =  0.  3  =  .5.  * 
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Figure  20:  Domain  for  oomputation  ot  hourviar 
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Figure  22:  One-dimensionai  tension  test. 
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(a)  20  x  4  elements. 


(b)  40  x  3  elements. 


(c)  80  x  3  elements. 
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Figure  27:  Restrained  panel.  Comparison  of  beam  theory  solution  with  finite  element  solii 
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Figure  32:  Solution  of  a  compression  corner  with  no  artificial  dissipation  and  lame  time  sten 

(CFL  =  10). 
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Figure  33:  Problem  statement  for  a  supersonic  boundary  layer  problem. 


Figure  34:  Temperature  profiles  for  a  supersonic  boundary  layer  problem 


Figure  35:  Problem  statement  for  a  tint  plate  pioblem. 


Figure  37:  Heating  rates  for  a  flat  plate  problem. 
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Figure  40:  Viscoplastic  stress  histories  in  the  aerodynamic  skin. 
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Figure  41:  Viscoplastic  stress  histories  in  the  heat  exchanger  fins 
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Figure  45:  Contours  of  a,j  stress  (a)  t  =  0.05s.  (b)  t  =  0.5s.  (c)  t 
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Aerodynamic  skin  with  fin  debond 


Damaged  fins 


(b)  Solution  obtained  with  large  displacement  theory. 


|  Figure  47:  Maximum  deformation  of  structure  ( t  =  0 .os)  with  damaged  heat  exchanger  fins 

Deformation  is  shown  to  scale. 
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Figure  4S:  Viscoplastic  deformation  of  aerodynamic  skin  for  undamaged  and  damaged  heat 
exchanger  fins. 


